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ABSTRACT 


A method is presented in satellite altimetry that attempts to 
simultaneously determine the geoid and sea surface topography with 
minimum wavelengths of about 500 km and to reduce the radial orbit 
errors caused by geopotential errors. 

The modeling of the radial orbit error is made using the linearized 
Lagrangian perturbation theory. Secular and second order effects are 
also included. After a rather extensive validation of the linearized 

equations, alternative expressions of the radial orbit error are derived. 
A Fourier series formulation allows for easier computations, examination 
of the frequency content of the error and computation of different 
statistics. A geographic representation with respect to geocentric 
latitude and longitude gives a significantly better insight on the spatial 
variations of the error. Similar expressions are derived for the geoid 
undulation error and sea surface topography. Numerical estimates for 
the radial orbit error and geoid undulation error are computed using 
the differences of two geopotential models as potential coefficient errors, 
for a Seasat orbit. 

To provide statistical estimates of the radial distances and the 
geoid, a covariance propagation is made based on the full geopotential 
covariance. Accuracy estimates for the Seasat orbits are given which 
agree quite well with already published results. 

Observation equations are developed using sea surface heights and 
crossover discrepancies as observables. A minimum variance solution 
with prior information provides estimates of parameters representing the 
sea surface topography and corrections to the gravity field that is used 
for the orbit generation. The potential of the method is demonstrated in 
a solution where simulated geopotential errors and the Levitus sea 
surface topography are used to generate the observables for a three 
day Seasat arc. The simulation results show that the method can be 
used to effectively reduce the radial orbit error and recover the sea 
surface topography. 
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CHAPTER I 


INTRODUCTION 


Until the middle of the last decade only about one third of the 
ocean surface was surveyed by sporadic oceanographic missions. With 
the advent of satellite altimetry a global observation system was 
established to provide synoptic observations of the oceans. The 
usefullness of satellite altimetry has been demonstrated by a series of 
altimeters of increasing precision flown on Skylab, Geos3, Seasat and 
Geosat. The Seasat altimeter in particular, has shown for the first 
time that spaceborne microwave altimeters are capable of obtaining 
global measurements of ocean topography with sufficient accuracy to 
be useful for oceanographic studies. Looking into the future, NASA 
and the French Centre National d’ Etudes Spatiales (CNES), in 
cooperation with other national partners, have been planning the 
coordinated Topex/Poseidon dedicated satellite mission for observing 
the oceans. 

The data collected from the different satellite altimeter missions 
have allowed for studies that led to substantial improvement in our 
knowledge of the gravity field of the earth, the shape of the oceanic 
geoid, ocean tides, current structure and oceanic behavior, crustal 
structure, solid earth dynamics and others. All these studies have 
produced a large volume and variety of scientific literature. One 
could mention, for example, three special issues of the Journal of 
Geophysical Research (Vol. 84, B8, 1979 on Geos3; Vols. 87, C5, 1982 
and 88, C3, 1983 on Seasat); Oceanography from Space, Marine Science 
(Vol. 13, 1981); Marine Geodesy (Vol. 8, 1984); and the Journal of 
Astronautical Sciences (Vol. 28, Oct. -Dec. 1980, on orbit determination 
for Seasat). Three comprehensive reviews of the work published in 
the USA alone from 1975 to date have been made by Stanley (1979), 
Marsh (1983), and Douglas et al. (1987), while a substantial number of 
publications is contained in the proceedings of a number of 
international symposia. 

From the very early days of Geos3 it has been clear that 
altimetric measurements of the oceans have many sources of error 
which are a complex function of position and scale. The most 
important of these errors were the ones associated with the orbit 
determination of the satellite and the altimeter itself. Continuous 
technological improvements have substantially improved the precision 
of the altimeter (10 cm noise for Seasat and even lower for future 
missions) so that altimeter noise is not considered as a serious error 
source anymore. The orbit determination related errors result from 
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the inability to accurately compute the radial component of the orbit 
at the time of the altimeter observation. This radial error is due 
primarily to errors in the tracking to the satellite, incomplete 
tracking, and errors in the modeling of the gravity field, air drag 
and solar radiation pressure. 

By far, the major component of the radial orbit error is due to 
the gravity field errors. Gravity field models and tracking of the 
70’s have resulted in radial errors of about 5 meters for Geos3. With 
the implementation of the GEM9 gravity field and use of altimeter 
observations, radial errors were reduced to the 2 meter level (Lerch 
et al., 1979). A similar improvement has been reported for the Seasat 
orbits. The first calculations of the radial component of the Seasat 
orbits had an error of about 5 meters which was reduced to about 1.5 
meters in the ephemeris that was used for the generation of the 
Geophysical Data Records. This error was further reduced to about 
70 cm in the early 1980’s (Lerch et al., 1982). 

The ephemeris accuracies for both Geos3 and Seasat provided 
accuracies of the sea surface that were prohibitive for many 
investigations, such as determination of sea surface topography and 
of an accurate oceanic geoid that can be used together with 
terrestrial data, for the computation of high degree and order gravity 
fields. So both Geos3 and Seasat observations needed to be analysed 
to reduce the radial error. New procedures were developed, ranging 
from fits of the altimetric surfaces to some long wavelength geoid 
computed from a satellite derived gravity model, to crossover analysis 
and optimal filtering. For the Seasat data analysis the most dominant 
method was the one that used crossover discrepancies as observables. 
The radial error, being recognized of long wavelength nature, was 
modeled either by a bias and tilt, or by a Fourier series of some low 
frequencies. These types of analyses have reduced the radial error 
to a reported level of 20-30 cm (e.g. Rowlands, 1981). Further 
combination with adjusted Geos3 data has produced accurate and 
detailed maps of the ocean surface (Marsh et al., 1984, Rapp, 1985, 
1986). Additionally the adjusted Seasat data with the combination of 
the GEML2 gravity field (Lerch et al., 1985a) have produced the first 
qualitative maps of the long wavelength permanent sea surface 
topography and ocean circulation (Tai and Wunsch, 1983, 1984, Douglas 
et al., 1984, Engelis, 1983,1985). 

Recently it has been recognized that in order to efficiently 
remove the radial errors from altimeter data, a more analytical 
modeling of these errors is necessary. Such a motivation was 
initiated by the early work of Anderle and Hoskins (1977) who showed 
that the errors are not random in nature but geographically 
correlated. Subsequent studies by Colombo (1984) and Wagner (1985) 
made use of the Lagrangian perturbation theory to derive expressions 
for the radial errors. These studies were followed by works of 
Tapley and Rosborough (1985), Rosborough (1986), Pisacane (1986), 
Mazzega (1986) and Tai and Fu (1986). Wunsch and Gaposchkin (1980), 
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Wagner (1986) and Marshall (1985) have proposed ways to 

simultaneously determine the permanent sea surface topography, 
improve the geoid and reduce the radial error. They have all used 
simulations to demonstrate the potential of their methods. 

The present study basically treats the same problem i.e., to 
analytically represent the orbit error due to the geopotential and 
describe its characteristics. Furthermore it attempts to optimally 
separate the long wavelength sea surface topography from the geoid 
at comparable wavelengths, while reducing the orbit error. This 
study has been independent and in parallel to most of the above 
mentioned investigations. A rather substantial use of Colombo’s 
results (Colombo, 1984) has been made and some of his notation is 
adopted. Some errors in his derivations have also been identified and 
corrected. This study extends the formulation to represent the geoid 
error and sea surface topography in a Lagrangian form. A Fourier 
series expression with respect to time and a geographic 
representation with respect to geocentric latitude and longitude are 
also developed for all the quantities. Then matrix expressions are 
developed and variances of the radial distances and the geoid as well 
as their error correlations are computed. The analytical formulation 
also allows for the development of observation equations to solve for 
parameters representing the gravity field and the sea surface 
topography. 

In a more detailed outline, this report is structured as follows. 
Chapter 2 makes a summary on the usefullness and applications of 
satellite altimetry as well as on solutions attempted in the past to 
remove the radial orbit errors. 

In Chapter 3 the linearized mathematical model to be used for the 
determination of the sea surface topography and the reduction of the 
orbit errors is developed. 

The detailed modeling of the radial orbit errors, geoid undulation 
errors and sea surface topography is made in Chapters 4, 5 and 6. 
Chapter 4 starts with the necessary background on the Lagrangian 
perturbation theory and the definition of a reference orbit to develop 
linear expressions for the radial error. An extensive validation to 
assess the accuracy of the theory is made by comparison with 
numerically integrated orbits. In Chapter 5 alternative equations of 
the radial error are developed. More specifically, a Fourier series 
expression with respect to time and a geographic representation with 
respect to geocentric latitude and longitude are computed, as well as 
expressions contained in Wagner (1985) and Rosborough (1986). 
Numerical "estimates" of the radial error and its properties for a 
Seasat orbit are then provided by assuming potential coefficient 
"errors" that are generated by differencing two satellite derived 
gravity fields. Chapter 6 analyzes the geoid undulation errors and 
sea surface topography along the lines of the previous Chapters and 
points out differences and similarities with the radial errors. 


4 


All these models are used in Chapter 7 to derive matrix 
expressions for all quantities and compute variances for radial 
distances and geoid undulations as well as their error correlations. 
This covariance propagation is made using the full covariance matrix 
of a gravity field. Numerical estimates on the accuracies of the radial 
distances for Seasat arcs are presented both along the subsatellite 
tracks and on a geographic grid. 

The simultaneous determination of the sea surface topography and 
the reduction of the radial errors using altimeter data is addressed in 
Chapter 8. Problems involved with the estimation are identified and 
crossover discrepancies, in addition to sea surface heights, as well as 
prior information for the parameters are used to improve the solution. 
A simulated solution for a three day Seasat arc is made to 
demonstrate the potential of the method if it is to be used with real 
altimeter data. 

In Chapter 9 conclusions from this investigation are drawn and 
recommendations for further work are given. 


CHAPTER II 


DESCRIPTION AND APPLICATIONS OF SATELLITE ALTIMETRY 


2.1 Introduction 

Satellite altimetry can provide the height of the sea surface 
relative to a reference ellipsoid. This height is a function of ocean 
variability and the geoid. Typically the geoid varies by up to a 
hundred meters relative to a reference ellipsoid with wavelengths 
ranging from a few to thousands of kilometers. Superimposed on the 
geoid are variations due to ocean dynamics. These effects are much 
smaller in amplitude than the geoid having magnitudes on the order 
of 1-2 meters. Both quantities are very important for geophysical 
and oceanographic studies. Therefore, before proceeding to the 
description and applications of satellite altimetry, a short discussion 
on geophysical and oceanographic uses of the geoid and ocean 
variability is appropriate. 

Knowledge of the oceanic geoid is very useful in geophysics, 
since one can infer information on seafloor topography, rigidity of the 
lithosphere and mantle convection. Changes in the shape of the geoid 
with maximum wavelengths of a few hundred kilometers are primarily 
caused by changes in the seafloor topography. Seamounts are 
particularly prominent, producing changes in the geoid from several 
centimeters to ten meters over distances of tens of kilometers. Thus, 
detailed and accurate maps of the geoid are needed to reveal the 
existence of seamounts. At longer wavelengths, variations of the 
geoid combined with other information (seismic and geochemical) can 
give information about the strength and thermal properties of the 
lithosphere. Finally the regional features of the geoid are 
attributable to deeper structures below the lithosphere and are 
related to mantle convection and other processes. 

Knowledge of the signature in the sea surface due to ocean 
dynamics and external forces is important because it can lead to the 
determination of the ocean circulation. Theory and observations have 
shown that to a first order of approximation, the large scale motion of 
the sea is geostrophic. Assuming a geostrophic balance and 

hydrostatic equilibrium, the equations of motion for the oceans at any 
depth can be developed. These equations can be solved using 

salinity and temperature observations from ships and considering a 
reference surface to be known. Such a reference surface can very 
well be the geoid, so its deviations from the sea surface (sea surface 
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topography) as well as the ocean surface velocities need to be known 
both spatially on a worldwide basis, and temporally over a long period 
of time (Coleman, 1981). 

The sea surface variations are composed of flows with different 
scales both spatially and temporally. The long wavelength spatial 
variations are basically time independent (stationary SST) and are 
attributed to the major oceanic currents. Time variations of SST are 
primarily of mesoscale spatial nature (eddies). Finally there are a 
number of smaller scale higher frequency phenomena such as internal 
waves, tides, tsunamis and others. These phenomena are not 
geostrophically balanced and are not part of the circulation itself. 

For several decades oceanographers have been trying to 
determine the sea surface topography and ocean circulation. Problems 
in these determinations are primarily due to the inability of 
determining a reliable reference surface and due to the shipboard 
data sampling, since existing observations are taken sporadically over 
different periods of time and so they do not represent a homogeneous 
behavior of the oceans over an extended period of Lime. The most 
recent estimates of ocean circulation and sea surface topography have 
been computed by Levitus (1982). 


2.2 Satellite Altimetry 

The principle of using an altimetric system for earth related 
studies is fairly simple. A radar microwave altimeter onboard a 
satellite measures its distance to the ocean surface. If the geocentric 
distance of the satellite is also known at the epoch of observation, 
then the height of the sea surface with respect to a geocentric 
reference ellipsoid is (Gopalapillai, 1974) 


h - 


P - - 


1 - 

l r 1 


i n~ 


' 2.11 


where p is the nadir distance from the center of mass of the satellite 
to the sea surface, h is the ellipsoidal height of the sea surface and 
r c , r are the geocentric distances of the ellipsoid and the satellite 
respectively. The last Lerm in equation (2.1) is a correction that 
accounts for the non colinearity of the geocentric distance r and the 
nadir distance p. It is a function of the eccentricity e of the 
reference ellipsoid and the geodetic latitude <t>,. of the satellite 
position. 

The nadir distance, p, from the center of mass of the satellite to 
the ocean surface is provided by the altimeter observation. Each 
altimeter observation is the outcome of considerable preprocessing of 
the raw measurement of the altimeter. The fundamental measurement 
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is obtained from the round trip travel time, At, required for an 
electromagnetic pulse, transmitted by the onboard electronics and 
emitted from the altimeter antenna, to reach the ocean surface, return 
back to the antenna and be detected and recorded by the onboard 
electronics. Sensor related corrections must be applied before the 
travel times are converted to distances (Hancock et al., 1980). These 
distances have further to be corrected for biases, reference to the 
center of mass of the satellite, relativistic effects and atmospheric 
propagation delays (Tapley et al., 1982). The latter can be separated 
into tropospheric and ionospheric effects. The retarding effects of 
the troposphere amount to about 2.5 meters with the dominant portion 
coming from the dry component. This effect can be modeled to the 
1-2 cm level, using surface pressure measurements. The maximum 
correction due to the wet tropospheric delays is 40 cm. These delays 
have rapid changes so the best way to correct for them is to include 
a scanning multichannel microwave radiometer (SMMR) on board the 
satellite that can measure the error caused by the pulse delays due 
to the water vapor. Ionospheric effects on the other hand can be 
computed to the 3 cm level by using a model compiled from 
independent data (Lorell et al., 1982) or by using a dual frequency 
altimeter. After these corrections are made, individual altimeter 
measurements can be averaged (over a couple of seconds) to produce 
one observation which is regarded to have an error being mostly 
white noise. Additional "in flight" calibration of the altimeter is 
needed (Kolenkiewicz and Martin, 1982; Marsh and Williamson, 1982) to 
assure that no bias exists in the altimeter observations. 

The Seasat instrumentation included a single frequency radar 
altimeter and a five frequency microwave radiometer. The reported 
accuracy of the observations was 10 cm for measurements averaged 
over a 1 sec interval. The footprint of each observation had an 
average width of 7 km along the satellite track on the surface of the 
earth. The Geosat altimeter is basically an improved version of the 
Seasat technology providing an accuracy of 5 cm after a 1 sec 
averaging. The Topex/Poseidon satellite is expected to carry 

instruments from both NASA and ONES. The NASA instruments include 
an advanced dual frequency radar altimeter and a three frequency 
microwave radiometer. After the necessary calibrations, the altimeter 
observations are expected to have an accuracy of 2 cm after 
averaging over a 3 sec interval which corresponds to an average 
footprint of 20 km along the satellite track. The ONES altimeter is an 
experimental single frequency solid state altimeter. It is expected to 
operate only 5 % of the satellite lifetime at regular intervals and 
provide observations with an accuracy of 4 cm after similar averaging 
(Stewart et al., 1986). 

The geocentric distance, r, of the satellite is determined from the 
computed ephemeris which is obtained by a high precision orbit 
determination process. The determination of the precise ephemeris 
involves the adjustment of the initial position of the satellite and 
several other parameters representing mainly coefficients of non 
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conservative forces. This adjustment minimizes in the least square 
sense the differences between the tracking data and the 

corresponding values that are computed using known station 
coordinates and comprehensive force models for the earth’s 
geopotential, lunisolar attraction, air drag, solar radiation pressure, 
earth tides, earth albedo, relativistic effects and others. The 
parameters to be adjusted are usually estimated at short intervals of 
a few days. Based on the initial positions and velocities for each 
interval, the orbits are then computed using a numerical integration 
of the force models. The solution is iterated until convergence. 

There are several errors involved in the determination of the 
precise ephemeris, associated with the observations to the satellite 
and the force models. Observation errors include errors in the 
tracking data, imperfect refraction corrections, errors in the station 
coordinates, and non global tracking data, both spatially and 
temporally. Tracking of altimeter satellites is usually done by laser 
ranging, Unified S-Band (USB) radiometry, and Doppler. For the 
Seasat satellite, laser tracking data had accuracies of 10 cm and 50 
cm for the Goddard Space Flight Center (GSFC) and the Smithsonian 
Astrophysical Observatory (SAO) lasers, respectively. The USB range 
rate data had a precision of 0.1 mm/ sec with an accuracy degraded 
by atmospheric effects. Finally the Doppler data, collected by the 
Department of Defense, had a precision of 10 cm at 30 sec integration 
intervals in range difference measurements. Spatial coverage was 

sparse for the laser data and almost global for the USB and Doppler 
data (Tapley and Born, 1980). 

The Topex/Poseidon satellite will have four modes of tracking. 
The NASA instruments of tracking, onboard the satellite, will be a 
laser retroref lector, a Tranet Beacon and a Global Positioning System 
experimental receiver. The CNES instrument is a radiometric tracking 
system using a one way Doppler system, called Doris, that receives 
signals transmitted from the ground. The primary tracking will be 
provided by the Defense Mapping Agency Tranet network and is 
expected to provide single pass radial accuracies of 13 cm. Laser 
tracking will be used primarily for the altimeter calibration and the 
satellite ephemeris validation. The observational accuracy of the laser 
ranging will be 2-5 cm. The Doris and GPS systems are expected to 
provide radial accuracies of 10 cm (Stewart et al., 1986). 

The most dominant errors in the precise ephemeris are due to 
gravity field modeling errors. Gravity field errors propagate into 
errors in the estimation of the initial state vector of the satellite and 
errors along the computed orbits from numerical integration. These 
errors have a variety of signatures and are systematic both spatially 
and temporally. As it will be seen later on, they are primarily of 
very long wavelength nature. Air drag and solar radiation pressure 
errors are smaller than gravity field errors. For altimetric satellites 
though they can be significant. The reason is that due to the 
complexity of their shape, the accurate modeling of the surface forces 
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acting on them is extremely difficult and so , only approximations are 
used in the modeling. Partial compensation for those approximations 
is made by adjusting, together with the initial state vector, the 
coefficients representing air drag and solar radiation pressure. 

Satellite derived gravity field models are being computed since 
the early 1960's and are continuously improved. Their improvement is 
reflected in the improvement of the accuracies in the computations of 
the radial distance to altimetric satellites. Based on prelaunch 

knowledge of the gravity field and other forces, the error in the first 
calculations of the radial component of the Seasat orbits was about 5 
m. After tailoring the gravity field with Seasat laser and USB 
tracking data and some Geos3 altimeter data, the error was reduced to 
about 1.5 m. This gravity field, PGSS3, was further improved with 
the incorporation of Seasat altimeter data to produce the PGSS4 
gravity field, which gave a radial orbit error of about 70 cm (Lerch 
et al., 1982). In order to meet the 10 cm requirements for radial 
position accuracies for the Topex/Poseidon mission, a new effort is 
currently underway at GSFC and the University of Texas at Austin to 
compute a very accurate gravity field. Preliminary solutions have led 
to the generation of the GEM-T1 field (Marsh et al., 1986b) which 
yields orbit accuracies much higher than any of the previous fields. 
Further improvement is underway since this model is expected to 
provide accuracies of 20-25 cm (ibid) for the Topex/Poseidon orbits. 

Considering the observational errors, A p, and the errors in the 
determination of the geocentric distance of the satellite, Ar, equation 
(2.1) becomes (after neglecting the correction term) 


h = r c “ P Q bs - r E + Ar - A/> (2.2) 

where p obs is the actual altimeter observation and r c is the computed 
geocentric distance. The sea surface heights, on the other hand, can 
be written as 


h = N+ e + f t + T + w 


(2.3) 


where N is the geoid undulation, f is the stationary SST, are the 
time variations of SST, r are the solid earth and ocean tides and w 
are other motions of the oceans like wind driven waves and others. 
With the implication of (2.2), (2.3) becomes 


r c “ Pobs - r E + Ar-Ap-N + f + f t + T + w 


(2.4) 


or, considering that h c = r c - p obB - r E 
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h c = N+ f + f^ + r + w-"Ar+ A p 


(2.5) 


where h c are the "observed" sea surface heights. Estimates of geoid 
undulations in (2.5) can be computed by using a set of potential 
coefficients up to a certain degree and order. Then 


h c ~ N„ + AN C + AN° +f+ft+T+w— Ar + A p 


( 2 . 6 ) 


where N 0 is the computed undulation, AN C is the undulation error due 
to errors in the modeled gravity field (commission error) and AN 0 is 
the undulation signal that iB not modeled from the given field 
(omission error). Setting 


Ah = h c - N 0 (2.7) 

equation (2.6) becomes 


Ah = AN C + AN° +f+f t +r+w-Ar+Ap 


( 2 . 8 ) 


where Ah are the so called residual sea surface heights that can be 
used in stationary SST determinations. 

Of all the quantities contained in the observed sea surface 
heights, the geoidal undulations have by far the largest magnitude 
with a maximum of about 100 meters and a global root mean square 
variation of about 30 meters. They are primarily of long wavelength 
nature (over 90% of the geoid signature is at wavelengths longer than 
1000 km). High frequency variations with wavelengths down to 10 - 
20 km and amplitudes ranging from several centimeters to tens of 
meters do exist (Rapp, 1985). The stationary SST, as shown by 
oceanographic evidence, is also of long wavelength nature with 
additional high frequency contributions in energetic regions of the 
oceans. The time variations of SST and the tides are mesoscale 

phenomena while the orbit errors are long wavelength. 

For oceanographic applications, mesoscale SST variations can be 
determined by differencing repeat altimeter tracks. Any long 
wavelength signatures in these differences are attributed to long 
wavelength orbit errors and can be easily removed (Cheney et al., 
1983). On the contrary, as shown from equation (2.8), the 

determination of the stationary SST requires the knowledge of an 
accurate geoid (very small AN C and AN°) and the effective reduction 
of the orbit errors. 
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For geophysical applications and particularly for studying the 
local features of the oceans (e.g. seamounts) the long wavelength 
small oceanographic signal in equation {2.5) can be treated as noise 
since it does not have any effect in the high frequency spectrum of 
the geoid. This is not always true for the radial error as it can be 
seen later on and so it must be removed. For the determination of an 
accurate geoid at longer wavelengths the separation of the geoid from 
the stationary SST and the reduction of radial error become critical 
since both quantities have signatures in comparable wavelengths. 

So the problem in satellite altimetry is a combined one requiring 
the simultaneous determination of both the geoid and stationary SST 
while reducing the radial error. The observables that can be used 
for this purpose are the sea surface heights which are linearly 
related to all the above quantities. 


2.3 Existing Solutions of the Combined Altimetric Problem 

Solutions to the combined altimetric problem that have been 
attempted in the past have focused on the reduction of the radial 
errors and the determination of a time averaged sea surface. Then, 
from that adjusted sea surface, together with the implementation of an 
independent satellite derived gravity model, estimates of the 
stationary SST were derived. 

The method that has been primarily used for the determination of 
the radial errors is the minimization in the least square sense, of the 
crossover discrepancies observed at the intersections of the 
ascending and descending groundtracks of the satellite. These 
discrepancies do not contain any geoid and stationary SST 

information. They only represent radial errors and SST time 

variations as well as the time variable part of mismodeling of 

quantities such as tides. In using the crossover discrepancies it was 
assumed that the radial error is independent in ascending and 

descending arcs. So minimization of the discrepancies would also 
effectively minimize the radial error. 

Based on the above principle for the solution, several techniques 
have been developed. In one technique followed by Rowlands (1981) 
for the Seasat sea surface heights, the Seasat arcs were broken into 
short arcs bounded by the ocean boundaries. For each of the short 
arcs, modeling of the radial errors was made by a bias and a tilt or 

only by a bias depending on the length of the arc. In order to 

remove the time variations of SST, repeating arcs from the repeat era 
of Seasat were averaged and used for the adjustment. In order to 

provide a reference for the solution an arc of 35 minutes of length in 

the Atlantic Ocean was held fixed. Then a global adjustment of the 
averaged arcs was made to reduce the crossover discrepancies from 
1.5 meters to 0.28 meters. The adjusted arcs were then held fixed to 
adjust the arcs of the non repeat era of the satellite in regional 
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solutions. 

Similar solutions, with small variations of the above technique, 
were made by Cloutier (1981), Marsh and Cheney (1982), and Marsh et 
al., (1982, 1986a). In Marsh et al. (1986a) the ocean surface was 
divided into small diamond shaped regions. Then minimization of the 
crossover discrepancies between the arcs within those regions was 
made to solve for the radial error that was modeled as a single bias 
for each arc. The resulting surfaces were then used to create a 
unique surface by adjusting the biases and tilts among these 
surfaces. All these solutions that use the crossover discrepancies as 
observables, although they have improved the estimates of the sea 
surface, are subject to several problems. The most important one as 
will be seen later on, is that the radial error has a component that is 
common to both ascending and descending arcs and cancels out 
during the generation of the crossover discrepancies. Therefore it 
cannot be removed during the crossover adjustment and so it exists 
in the adjusted sea surface. Additionally all the error along the arcs 
that are held fixed propagates to the adjusted surface. Finally all 
the effects like mismodeling of tides and SST time variations (in the 
non repeat arcs) alias the solution. 

An alternative technique for removal of the radial error is the 
one that uses the residual sea surface heights as observables (see 
equation 2.8). The radial error, being recognized of long wavelength 
nature, is modeled by a Fourier series of low frequencies. Then, a 
solution for the Fourier coefficients yields corrections to the radial 
distance and eventually to the ocean surface. The potential of such a 
method to improve the satellite ephemeris was demonstrated by Goad 
et al., (1980) who also recognized limitations of the method because of 
aliasing due to the long wavelength geoid undulation commission 
errors and the existence of the oceanic effects. A similar solution 
was also used by Rapp (1979) who used the residual sea surface 
heights, in combination with crossover discrepancies to adjust the 
Geos3 data. 

Since the Seasat ground tracks are not dense enough to allow for 
the detailed mapping of the sea surface, as is needed for geophysical 
applications, the combination of Geos3 and Seasat sea surface heights 
was made, using basically the same procedures and considering all the 
crossovers created between the ascending and descending arcs of the 
two satellites (Liang, 1983, Marsh et al., 1984). These combinations 
created high resolution mapB of the sea surface and the gravity field 
of the oceans (Rapp, 1985, 1986) but they also introduced additional 
problems such as across track high frequency residual errors that 
can be easily detected in the maps. These residual errors are due to 
the fact that the non modeled part of the orbit error has different 
magnitudes in Geos3 and Seasat, resulting in a signature of long 
wavelength nature along track but of a high frequency nature across 
track. 
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The adjusted Seasat sea surface heights have been used to 
determine the stationary SST. In order to model the geoid recent 
satellite derived gravity models were used and residual sea surface 
heights were computed. These models have accuracies that only allow 
for the determination of very long wavelength stationary SST. Using 
the GEML2 gravity field, estimates of the stationary SST with 
wavelengths on the order of 3000 km and longer have been computed 
by Engelis (1985), Tai and Wunsch (1983, 1984) and Douglas et al., 
(1984) that agree quite well with oceanographic estimates at least 
qualitatively. 


CHAPTER III 


FORMULATION OF THE SOLUTION 


3.1 Introduction 

To overcome the problems that arise from the empirical solutions 
of the combined altimeter problem, a more general solution has to be 
made, involving the simultaneous determination of the geoid, 
stationary SST and other ocean effects and the reduction of radial 
errors. In order to do that, a mathematical model relating all 
quantities has to be developed, in which the dynamic properties of 
the orbit are to be considered. The basic equation that relates all 
quantities of interest is obtained from (2.1) and (2.3) 


h=r-p-r E =N+f+f t +T+w 


(3.1) 


In the context of this study which concentrates on the determination 
of the stationary SST, the sea surface time variations and the tides 
will be omitted from equation (3.1). This could in principle be a 
justifiable omission since sea surface variations can be determined 
independently (Cheney et al., 1983). Furthermore, the major ocean 
tides are obtained by contemporary ocean tidal charts, which are 
believed to be reasonably good (Schwiderski, 1980) although direct 
evidence for this is limited to data from tidal stations scattered 
widely along the coasts and islands. Recently a new procedure has 
been developed to also determine the minor tidal constituents by 
using a linear interpolation on the admittance function of the major 
ocean tides (Christodoulidis et al., 1986). Finally solid earth tides are 
considered to be known much better than ocean tides at present. So 
(3.1) reduces to 


h=r-p-r E =N+C (3.2) 


3.2 Mathematical Model 

It is well known that the state of a satellite at any epoch t can 
be represented by a six dimensional vector s, where three of the 
components of s represent the position and the other three the 
velocity of the satellite. Then s x is the initial state vector 
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corresponding to the beginning of the orbit which usually coincides 
with the time origin. Let p be an n-dimensional vector of parameters 
associated with the various forces acting on the satellite and q be an 
m-dimensional vector representing all other non force model 
parameters. Then any observation, f, to the satellite is generally a 
nonlinear function of those parameters 

f = f(s(p), q) (3.3) 


where f has a functional dependence on the parameters p through the 
state vector s. The task in orbit determination is to solve equation 
(3.3) for the state vector and/or the parameters. Since this is a 
nonlinear problem, initial values p D , q Q for the parameters as well as 
an approximate state vector s 0 are required for its linearization. 
Then a Taylor series expansion gives 


f = 


f " + H 45 


3f 

3q 


Aq 


(3.4) 


where f Q is the computed observation from (3.3) using the given 
approximate values, and As, Ap, Aq are small corrections to the state 
vector and the parameters. The corrections As and Ap are 
interrelated through the functional relationship s(p). As a matter of 
fact the state s of the satellite has a nonlinear dependence on the 
force parameters p as well as on the initial state s x as follows 

s = §(§!, p) (3.5) 


and so small variations As x and Ap result in small errors of the state 
vector 


. - 3s , _ ,3s . - 

As = As t + -TT7 Ap 

3s x 1 3p 


(3.6) 


where the partials represent matrices with elements the derivatives of 
the state at any time with respect to all components of s x and p 
respectively. Substituting (3.6) into (3.4) we obtain 


* * A 3s t 3f 3s ^ 3f A _ 

f f ° + 3f 317 1 35 AP 3q 


(3.7) 
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Equation (3.7) is the general linearized mathematical model that 
relates any satellite observation with the initial state vector, the force 
field and the other parameters. The partials 3f/3s and 3f/3q are to 
be computed on the approximate orbit s c using the initial values p 
q„ of the parameters. The partials 3i/3s x are the solutions of the 
variational equations and they compose the so called state transition 
matrix. Equation (3.7) can be solved in the least squares sense to 
determine the initial state vector and any parameters of interest. The 
adjusted parameters are then used to integrate the equations of 
motion and determine the satellite state vector at any epoch. The 
solution is iterated until convergence. 

Equation (3.7) can be used to set up the mathematical model for 
the solution of the combined altimeter problem. The basic observable 
in this case is the distance p between the sea surface and the 
satellite. From (3.2) we have 


p r - h - r E '3.8) 

Neglecting for the moment the; dependence of h on the parameters p 
and using (3.7) for p , we obtain 


3r 

3s 




3s 

3S X 


~’T 


3r 3s 
3s 3p P 


3r - 


h - r « 


(3.9) 


wWe = ITT 

has been used. Since the; radial orbit error in altimetric satellites is 
ou the order of a few meters and is primarily of long wavelength 
nature, equation (3.9) is a good linear api>roximation of the 
corresponding non linear equation of the type (3.3) when the initial 
values r 0 , s u , p 0 , q 0 are the ones that have been computed (r 0 , s,,), 
or used (p Q , q D ) in a precise orbit determination process. A solution 
of equation (3.9) can then be considered as a second step to a two 
step solution of the combined altimeter problem. In the first step, 
the satellite observations (laser ranges, USB range rates and Doppler 
range differences), a very accurate gravity field, surface force- 
models, a tidal model and a set of station coordinates are used to 
provide estimates of the radial distance r, the initial state vector s r 
and surface force parameters p. In the second step the altimeter 
observations are used, together with r and s x , to improve the gravity 
field and the surface force parameters and to provide a more accurate 
sea surface. 
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For the purpose of this study we neglect the parameters 
representing the surface forces. Then the force parameters p 
represent only the potential coefficients of the gravity field. 
Considering also that there is no explicit relationship of the 
observables with the parameters q, we obtain 


P ~ 


r„ 


3r 3s - dr 3s 

3s 3sj S * 3s 3p 


4p - 


h - 


r E 


(3.10) 


The implementation of the stationary SST and geoid undulations in 
equation (3.10) can be made by using (3.2). The use of a reference 
gravity field for the geoid undulations gives 


h = N 0 



+ AN 0 + 


K. 

3p T 


Ap T 


r E 


(3.11) 


where p' is an n-dimensional vector containing potential coefficients 
of the reference gravity field, AN° are the omission errors and p T is 
a k-dimensional vector containing parameters of the stationary SST. 
For the SST representation a spherical harmonic expansion is also 
adopted although it has not been established yet that this is the best 
representation. A reference field could in principle be used for the 
stationary SST, using oceanographically derived estimates. But since 
its magnitude is on the order of 1-2 meters no such reference is 
necessary. Combination of (3.10) and (3.11) gives 


XT . 3r 3s . - 
P r o N 0 + a - a - As, 


3r 3s - 3N_ 

35 3p AP 3p' 


A P' 


" 1^7 A P T “ AN ° " r E 


(3.12) 


A very important aspect of equation (3.12) is the choice of the 
reference gravity field for the geoid undulations. The efficient 
mapping of the geoid requires a high degree and order gravity field 
model. This is in contrast with the satellite perturbations which are 
only sensitive, due to satellite altitude, to gravitational variations that 
can be described by a relatively low degree field. In principle one 
can use an existing high degree reference field (e.g., Rapp and Cruz, 
1986) so that the omission errors are small and solve for corrections 
to that field. This choice though can be impractical computationally 
and certainly will not provide good high degree estimates since 
altimetric observations are only taken over the oceans. So unless 
some combination with terrestrial gravity data is made, any attempted 
solution will seriously suffer from leakage effects. On the other hand 
high degree solutions can be obtained by techniques like satellite to 
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satellite tracking which are global in nature. So for the solution of 
(3.12), the reference field that is adopted for the undulations can be 
the same as the one used for the orbit determination. In such a case 
the omission error can be quite large and some special filtering has 
to be made to the residual sea surface heights Ah to remove it. This 
filtering will also remove all the stationary SST signature with 
harmonic degrees greater than the maximum degree of the reference 
gravity field. More details on the geoid representation, the omission 
errors and their removal are contained in Chapter 6. Assuming that 
AN 0 can be effectively removed, (3.12) can be written as 


P - P o 


3r 3s 
3s 3ij 


Asj+ 


(|r||_|N] A - _|£_ 

Ids dp dpJ dp T 


Ap T 


(3.13) 


where Po = r 0 - N 0 - r E 


(3.14) 


In terms of the residual sea surface heights equation (3.13) is 
written as 


Ah + 


3r 3s 
3s 3gj 


As j+ 


(3r 3s 3N1 . at . 

Idl 3£ “ 3^J 4p - afc Ap T = 0 


3p T 


(3.15) 


Equations (3.13) and (3.15) are the basic mathematical models of 
the combined altimeter problem. In order to be able to solve for the 
corrections As^ Ap and Ap T , the functional representation of the 
radial orbit error, undulation error and stationary SST has to be 
developed. This modeling will give insight into the properties of the 
above quantities and their correlations and will provide the 
observation equations for the solution. After the corrections are 
estimated the same models can be used to compute the radial orbit 
error, undulation error and stationary SST to provide a solution for 
the combined altimeter problem. These models will be developed in 
detail in Chapter 4, 5, and 6. 



CHAPTER IV 


ANALYSIS OF RADIAL ORBIT ERRORS 


4.1 The Satellite Orbit 

The state of a satellite at any epoch can be defined in a 
Cartesian coordinate system with origin at the center of mass of the 
earth. The x axis of such a system points towards the vernal equinox 
and the z axis is directed towards the celestial pole at some 
fundamental mean epoch. The y axis is perpendicular to the xz plane 
so as to form a right handed orthogonal system. This Cartesian 
coordinate system is semi-inertial since it is accelerated, together with 
the earth, from external gravitational attractions of the sun, the moon 
and other planets. In practice the orientation of this system is 
defined from quasi-stellar radio sources while its origin is defined by 
satellite orbit dynamics. An alternative but equivalent system can be 
defined if agreed upon values of precession and nutation are adopted. 
In such a system the z axis coincides with the spin axis of the earth 
that is consistent with the adopted nutation theory, the x axis is 
directed towards the vernal equinox and the xy plane is 
perpendicular to the z axis and defines the equator of the earth. 
The Cartesian formulation does not lend itself to any singularites in 
describing the state of a satellite and is very convenient in 
developing the equations of motion and integrating them. Alternative 
formulations also exist that are better suited to describe properties of 
the orbit. The most common one, that makes, the description of a 
satellite orbit particularly simple, is the formulation in Keplerian 
elements. The Keplerian elements define an ellipse that has one focus 
at the geocenter and is always tangent to the true orbit. That is 
why the ellipse is called an instantaneous or osculating ellipse and 
the corresponding elements instantaneous or osculating elements. 

The size and shape of the ellipse are defined by the semimajor 
axis, a, and the eccentricity, e, while its orientation with respect to 
the Cartesian system is defined by the inclination, i, the right 
ascension of the ascending node, 0, and the argument of perigee, w. 
The position of the satellite on the osculating ellipse is specified by 
the true anomaly, f, or equivalently by the eccentric anomaly, E, that 
is related to the mean anomaly, M, by Kepler’s equation 


M = E - esinE 


( 4 . 1 ) 
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which is a transcendental equation in E and can be solved iteratively- 
given M and e. The six Keplerian elements are equivalent to the 
rectangular components of position and velocity. They have the 

advantage, with the exception of the mean anomaly, that their 
variation in time is much smaller than the variation of the rectangular 
coordinates. The rate of change of the mean anomaly on the 
osculating ellipse can be approximated by the mean motion, n, which 
according to Kepler’s third law, is 


n = ^ M(t) = (pa" 3 )* (4.2) 


where p is the product of the universal constant of gravitation, G, 
times the mass of the earth and the atmosphere, M e . 

The motion of the satellite can be traced on the surface of the 

earth by the so called groundtrack of the satellite which is the line 

connecting all the subsatellite points. The satellite orbit can be 
tuned so as to form a periodically repeating groundtrack. This type 
of orbit is called a frozen orbit, and is used in altimetric missions. 
Altimeter satellites are placed in orbits of small eccentricity to keep 
the distance to the ocean surface always close to the optimum range 
of the onboard instruments. They are also put high to reduce air 

drag, but sufficiently low to have short orbital periods and finely 

spaced groundtrack. Their inclination is chosen so that most of the 
oceans are sampled. These requirements result in orbits with an 
altitude of about 1000 km, a small eccentricity (e ~ 0.001) and 
inclinations larger than 60*. The orbital periods are all close to 100 
minutes or about 14 revolutions per day. 

The relationships between the rectangular coordinates, the 
Keplerian elements and the geographic coordinates of the groundtrack 
will be used later on in this study and are given in Appendix A. A 
detailed description of the orbital geometry can be found in Kaula 
(1966). 


4.2 Equations of Motion 

The Newtonian equations of motion for the center of mass of a 
satellite in an inertial coordinate system can be written 


r = P 


(4.3) 


which is a set of three second order differential equations of the 
position vector r. P is the total acceleration vector and contains 
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accelerations due to the gravitational field of the earth, third body 
attractions, aerodynamic forces, solar radiation pressure, thrusting of 
the spacecraft rockets, attitude control system corrections and others_. 
If we consider only the gravitational acceleration of the earth then P 
is the gradient of the gravitational potential and (4.2) becomes 


r = vv 


(4.4) 


In order to express the equations of motion from rectangular 
coordinates to Keplerian elements the three second order equations 
have to be transformed to six first order equations. This 
transformation is made by Lagrange and is described in Kaula (1966). 
The Lagrangian equations of motion are 


da _2 3F 

dt na 3M 


de _ 1-e 2 3F _ (1-e 2 )* 3F 
dt na 2 e 3M na 2 e 3« 


di _ cosi 3F _ 1 3F 

dt na 2 (l-e 2 )*sini 3« na 2 (l-e 2 )*sini 3Q 


do cosi 3F (1-e 2 )* 3F 

dt na 2 (l-e 2 )^sini 3i na 2 e 3e 


dQ _ 1 3F 

dt na 2 (l-e 2 ) ,i sini 3i 

dM _ _ 1-e 2 3F 2 3F 

dt na 2 e 3e na 3a 


where the forcing function F is 



(4.6) 


The function R is known as the disturbing function and contains all 
the terms of the gravitational potential V except the zero degree term. 
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Inserting (4.6) in (4.5) we obtain the same form of equations, with F 
replaced by R, except for the last equation that becomes 


dM _ _ l-e a 3R _ _2 3R 

dt n na 2 e 3e na 3a 


Equations (4.5) and (4.7) have the general form 


s = s(s, p) 


(4.8) 


where the dependency on s is straightforward while the dependency 
on p is given through the disturbing function R. Obviously small 
variations in p result in small errors in the rates of the Keplerian 
elements 


. - 3s . - 
As = Ap 
3p 


(4.9) 


4.3 The disturbing function 

The disturbing function R which is equal to the gravitational 
potential V without the zero degree term, is a time invariant function 
in an earth-fixed geocentric coordinate system. It has the very well 
known form (Heiskanen and Moritz, 1967) 




£ £ _ 

X (C. cosmX + S. sinmX)P. (sin«t>) 

o *(R % m 


(4.10) 


where r,4»,X are geocentric distance, latitude and longitude in an earth 
fixed equatorial system, a e is the semimajor axis of the e_arth, Ci m and 
Sj. are the fully normalized potential coefficients and Pi m (sin4>) are 
the fully normalized associated Legendre functions. A more compact 
expression for the disturbing function can be obtained by considering 
the surface spherical harmonics which can be written in the form 


Y. (♦.X) 

*ma 


I P t (sin4>)cos(mX - ?a) 

a — o ^ 


(4.11) 
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Then (4.7) becomes 


B “® iU®;) j. .i. 5 i.. 


x) 


(4.12) 


where C. = C. and C. = S. 

*mO #mi *m 


The disturbing function R, in order to be used in equations (4.5), 
has to be transformed from the earth fixed coordinate system to an 
inertial coordinate system. This is a transformation that has been 
carried out in detail by Kaula (1966). The expression for the 
disturbing function R in terms of osculating Keplerian elements is 


R = H . 
a t -. 


<o ( - \ i t i to 

i ra i z f. (i) i 

= 2 k m=0 p ~ O * m P q— — o 


m-o p 


G. (e) S. («,M,0 ,0) 

*pq *mpq 


(4.13) 


where 


tmpq 


-s, 


fs u 


4 — m even 

cos[(4-2p+q)M + (4-2p)« + m(Q-0)] 

t— m odd 


1 4-m even 

sin[ (4-2p+q)M +(4-2p)u + m(Q-0)] 

4—m odd 


(4.14) 


and Fi rop (i) are the inclination functions, Gf pq (e) are the eccentricity 
functions and 6 is the Greenwich sidereal angle. A more compact form 
given by Colombo (1984) is 


R = 5 J Z X N 1 C. I Z F. (i) G (e) C. ( W ,M,n,8) 

a 2 m = 0 a =o l aJ p — o q-— « *mp *pq *mpqa 


(4.15) 


where 


C. = cos { ( i-2p+q)M + (4-2p)« + m(O-0) - 7 j[ a + ±( 1- ( -1 ))*-"] } 

«mpqa 4 4 


tmpq 


(4.16) 
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Equations (4.13) through (4.16) are expressed with respect to 
unnormalized coefficients. Similar expressions can be obtained if the 
fully normalized potential coefficients and the normalized inclination 
functions are used. The summation limit for the harmonic degree £, 
being infinity in theory, is in practice limited to a rather low value. 
Indeed, due to the altitude attenuation factor (a e /a) £ the higher the 
harmonic degree of a spherical harmonic coefficient, the lesser its 
effect on the potential and the gravitational acceleration acting on the 
satellite. It turns out that for a satellite with an altitude of about 
1000 km the series for R can be truncated at * max a 40. 


4.3.1 Inclination and Eccentricity functions 

The inclination functions Fi mp (i) are the result of a transformation 
of the Legendre functions from a plane parallel to the equator to a 
plane coinciding with the orbit plane. The expression for computing 
the inclination functions is (Kaula, 1966) 


(i) = I (21 2t)! (sini)'*~ m ~ 2t 

t t ! ( jC— t) ! (i-m-2t) !2 2i!-2t '• sini; 


&mp 



cos s i X 

C 


r*; 2t+a l (£! c ) <-» 


k 


(4.17) 


where k is the integer part of (i-m)/2, t is summed from 0 to the 
lesser of p or k, and c is summed over all values making the binomial 
coefficients non zero. 

This series representation of the inclination functions makes their 
computation quite costly or even impossible, due to hardware 
limitations, when the maximum harmonic degree is high. An efficient 
method to compute the inclination functions for practically any 
maximum harmonic degree has been developed by Goad (1987). This 
method uses the fact that the inclination functions (and their 
derivatives) are the Fourier coefficients of a series of values of 
surface spherical harmonics along a great circle of inclination i, and 
therefore uses a Fast Fourier Transform algorithm to recover their 
values. The method is stable, does not rely on short or long-period 
characteristics of the surface harmonics and does not require large 
amounts of computer resources even at degrees 180 and higher. An 
alternative technique that uses recurence relations has been 
developed by Kostelecky (1985). The technique is very efficient in 
computation time but requires large amounts of computer memory. 

The eccentricity functions, also known as Hansen coefficients, are 
needed to substitute r and f by a, M and e to obtain the final 
expression (4.12) for the disturbing function. Expressions for the 
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eccentricity functions are very complicated and can be found in Kaula 
(1966). For small eccentricities these functions become proportional to 
e*i' so for nearly circular orbits (e < 0.005) one usually retains only 
terms of q = -1, 0, 1 in the summation of equation (4.12). For small 
eccentricities the following approximations are valid 


°i,. = 1 + °< e > 



V, = 0 + °<*> 


a: 


P 1 


31-4P+1 

2 


(4.18) 


a -l+4p+l 
■e p -x “ 2 


where G' denotes the derivative of G with respect to e. The efficient 
computation of the eccentricity functions and their derivatives is 
being implemented by Goad (1987) in a method similar to the one for 
the inclination functions. 


4.4 The Radial Orbit Error 

From equation (3.13) the radial orbit error Ar can be expressed 
as 


. 3r 3s . - 
Ar = -r= -r=- ASt 
3s 3sj 1 


, 3 s 
+ TT r: Ap 
3s 3p 


(4.19) 


where the partial derivatives of r and s are to be taken along the 
precise orbit that is computed during the first step of the solution. 
Equation (4.19) can also be written as 


Ar = Ar x + Ar G 


(4.20) 


where 
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4rc = |f || ip (4.21) 


4r > = air 45 * (4 - 22) 

are the errors of direct gravitational origin and errors due to initial 
state vector errors. The explicit form of the radial orbit errors, of 
either origin, with respect to Keplerian element errors can be 
obtained by considering that the radial distance is 


r = a(l-ecosE) (4.23) 

Expanding E as a function of M, we get, to the order of e 3 
(Smart, 1977) 

cosE = (1— ^ e 2 )cosM - ^ ecos2M + ^ e 2 cos3M (4.24) 


So 


r = a(l-ecosM) + 0(e 2 ) 


(4.25) 


Then an error in a, e, M of either gravitational origin or due to the 
initial state vector will result in an error in the radial distance 


Ar = t-=- As = Aa(l-ecosM) - aAecosM + aeAMsinM 
3s 


(4.26) 


or 


Ar = Aa - (aAe+eAa)cosM + aeAMsinM 


(4.27) 


In the following, the radial orbit error from both initial state and 
gravitational origin will be analyzed. 
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4.5 The Radial Orbit Error of Gravitational Origin 


4.5.1 Linearized Equations of Motion 

The equations of motion of the type (4.8) or equivalently of the 
type (4.4) can be numerically integrated with a suitable integration 
technique. A description of several numerical integration techniques 
that are used for this purpose is contained in Cappellari et al., 
(1976). Equations (4.8) can also be integrated analytically by 

considering the following linearization 


s = §(§„) + || J s(i Q ) dt 


(4.28) 


where I 0 is a vector of approximate elements defining a reference 
orbit. The first term in equation (4.28) can be analytically integrated 
using the linearized Lagrangian perturbation theory. Let §i be this 
first order approximation to the osculating elements 


Si = / s(s 0 ) dt 


(4.29) 


Then the second term of (4.28) can, in principle, be also analytically 
integrated to provide the second order terms representing the 
coupling effects between the Keplerian elements. This analytical 
integration technique with different variations regarding the 
completeness in the representation of the disturbing function and the 
choice of variables (Keplerian, Delaunay, etc.) has been used in the 
past by Kozai (1959), Brouwer (1959), Brouwer and Clemence (1961), 
Kaula (1966) and others. But even with the incorporation of second 
order coupling effects, the accuracies with which the satellite state 
vector was computed were rather limited and so the analytical method 
has been replaced by numerical integration techniques. 

For the solution of the combined altimeter problem the quantities 
of interest are the small potential coefficient corrections giving rise 
to radial orbit errors on the order of 1-2 meters. For such 

magnitudes the linearized theory is expected to be quite accurate and 
is the one that is applied in this study. Consistent with the 

linearization of (4.8), (4.9) becomes 


ai *- 
As = ~r=- Ap 
3p 


6 = - 3 . s( § n ) 

3p 


_ 3s 3ij 
lp + 31 3f 


Ap 


(4.30) 
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or 


♦ • 0 g 

As = As(s 0 ) + •£! As! 


(4.31) 


where 



(4.32) 


is the first approximation to the Keplerian element errors computed on 
the reference orbit s 0 . Then (4.31) is integrated to 


As = As, 



Aijdt 


(4.33) 


Substituting into (4.21) we obtain 


. 3r *- 3r f 3s _ 

Ar G = af A8 i + gf J jf As, 


dt 


(4.34) 


Note that the partials dr/35 still need to be taken on the precise 
orbit. The first term in (4.34) represents the first order 
approximation of Ar G while the second term accounts for second order 
effects which arise from the interactions of the first order errors 
with all the potential coefficient terms. Before we proceed to the 
determination of the detailed expressions for the first and second 
order terms of the radial orbit error, the reference orbit s„ has to be 
defined. 


4.5.2 The Reference Orbit 

The major force acting on the satellite is the one due to the 
central force field of the earth. In such a field the satellite follows 
an elliptical orbit that obeys the three laws of Kepler. The size, 
shape and orientation of this ellipse are fixed and described by 
Keplerian elements that have no variations in time. Because the 
earth’s potential is not represented by only a zero degree term 
additional forces acting on the satellite make it depart from its normal 
orbit. The largest of those forces is due to the equatorial bulge of 
the earth which can be described by the (2,0) term of the 
gravitational field. The time independent part of the equations of 
motion, considering only the central force and equatorial bulge of the 
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earth are given by Kaula (1966) 


da 

dt 


de 

dt 


di 

dt 


0 

0 

0 


dcj 

dt 


dO 

dt 


dM 

dt 


3nJ 2 ag 
4(l-e 2 ) 2 a 2 

3nJ 2 ag 

2(l-e 2 ) 2 a 2 


n + 


3nJ 2 a? 


uiiu 2 q p 

4(l-e 2 ) 3 / 


[l-5cos 2 i] 


cosi 


JZI (3cos 2 i-l) 

d 


(4.35) 


where J a is the dynamic form factor of the earth. Then the resulting 
orbit will have the following elements 


a Q (t) = a 0 
e 0 (t) = e 0 
i Q (t) ~ i 0 

(4.36) 

«o(t) = u 0 + w(t-t 0 ) 

o 0 (t) = fi 0 + 6(t-t 0 ) 

M 0 (t) = M 0 + M(t-t 0 ) 

where (a 0 , e 0 , i 0 , « 0 , Q 0 , M 0 ) are the time invariant or mean Keplerian 
elements and t 0 is the time at the beginning of the orbit. Equations 
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(4.36) define an ellipse that is gradually changing with a secular 
motion. Its variations are a steady slow turning of the orbit plane 
due to a secular change in Q (precession of the line of nodes), a 
rotation of the major axis in the orbit plane about the geocenter due 
to a secular change in « (precession of the argument of perigee) and 
a slight departure from Kepler's third law given by equation (4.2). 

In order to use equations (4.36) to define the reference orbit, the 
mean Keplerian elements need to be determined. One method that has 
been used by Brouwer (1959) and Kozai (1966) is the following: since 

the largest periodic perturbations are the ones due to J 2 , a good 
approximation to the mean elements can be obtained if the J 2 periodic 
perturbations are computed and subtracted from the initial elements 
5j that are given from the orbit adjustment. These perturbations are 
computed in Appendix B. Then the mean values of the semimajor axis 
and the eccentricity, for example, are 


3 a 2 

a o = a - -x 3 2 ~ B ~ sin 2 icos2(o+M) (4.37) 

£ SI 

e 0 = e - J 2 [^] sin 2 icos(3M+2«) + ^ sin 2 icos(M+2«) (4.38) 


+ 


(| - | sin 2 i] cosm] 


where the elements on the right hand side are the initial elements. 

A more accurate determination of the mean elements can be 
obtained when a precise numerically integrated orbit is available for a 
time period that is an integer multiple of the periods of all orbital 
perturbations. The mean elements (a 0 , e„, i„) and the mean rates (M, 
«, fi) can be determined by averaging the computed osculating 
elements and rates over the time span for which the orbit is 
computed. Then (M 0 , « 0 , fl 0 ) are the mean values such that they give 
the best fits to the corresponding osculating elements in the sense s Q 
+ s(t-t 0 ) ~ s. This approach of determining the mean elements has 
the advantage, over the previous one, that in addition to the periodic 
effects due to J 2 , it also removes all the periodic effects due to the 
other terms of the geopotential. Additionally the secular rates contain 
all the secular (from all even zonal terms and not only from J 2 ) and 
very long period effects. In reality the two procedures are almost 
identical since the perturbations, other than from J 2 , are quite small 
and so the resulting mean elements are almost the same. This will be 
seen in Section 4.7 where a numerical validation of the theory 
developed in the present Chapter will be made. 
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For altimetric orbits, which are almost circular, the determination 
of the mean elements becomes more complicated. For such orbits the 
position of the perigee cannot be easily defined. As a result, 
variations in both the argument of perigee and mean anomaly are 
markedly non linear, while the eccentricity variations are not 
sinusoidal anymore. The perigee does not precess but oscillates about 
a mean value of *90* (Cook, 1966). It turns out that the mean 
anomaly has similar variations, so the total variations of (w+M) behave 
well and can be easily determined. The mean elements e and « can be 
computed using the theory developed by Cook who used the following 
non singular transformation 


f = ecos u 
7) = esincj 


(4.39) 


Cook has determined perturbations in f and 7j using the zonal only 
field and concluded that there are no perturbations of zonal origin in 
e and u when these elements are equal to 


7T 

w 0 “ 2 

(4.40) 

c 

e ° = 5 

(4.41) 

where 



c = 


^max 

n I 

£-3 


* w 


(J-D 

i(*+l) 


Pii(o)Pii(cosi) 


(4.42) 


k = — ^ nJ 2 (l-5cos 2 i) 


(4.43) 


i max being the maximum harmonic degree of the summation which is 
only over odd degrees. Pii(cosi) are the associated Legendre 
functions given by 


Pii (cosi) 


sini v (2l-2t) ! (~l) t 
~2 ir t (i-l-2t)!(l-t)!t! 


(cosi) 1 ® 1 2t 


(4.44) 


32 


the summation over t being from 0 to (i-l)/2. 

This type of orbit with u 0 = tt/ 2 and e 0 given by (4.41), (4.42) 
and (4.43) is the orbit that is desirable in satellite altimetry since it 
has the property of forming a very precisely repeating groundtrack. 
In reality « is very close to zero but not equal to it. For small 
values of u the mean values of e and u can still be computed using 
Cook’s theory. 

In the particular application of the combined altimeter problem the 
following reference orbit will be used. The mean values i ( ; > « 0 , Q 0 , M„ 
are identical to the initial values; a 0 , e 0 are computed by equations 
(4.37) and (4.41); M and Q are computed by (4.35) and « is computed 
by (4.35) supplemented by the effect of all odd zonal harmonics 
(consistent with Cook’s theory). It turns out (see Section 4.7) that 
such a reference orbit is very close to the one computed by 
averaging the osculating elements of a precise ephemeris. 


4.5.3 First Order Radial Orbit Error 

From (4.34) the first order radial orbit error is 


. 3r 

Arex = af ASl 


and written explicitly 


(4.45) 


Ar G1 = Aa x - (aAe l +eAa 1 )cosM + aeAM x sinM (4.46) 

where a, e, M are the osculating elements given by the precise 
ephemeris and Aa l( Ae x , AM X are computed analytically using the 
linearized perturbation theory. Expression (4.46) is valid in the 
general case of an elliptic orbit. For nearly circular orbits though 
like the altimetric ones, Ae x and AMi cannot be computed unless the 
non singular transformation (4.39) (applied for the errors) is used to 
compute Ae x and A« x . Then by computing A <a x + AM X the errors AM X 
can also be obtained. This transformation can be avoided by 
introducing the variable 


u = ecosM (4.47) 

which is also a non singular variable and can be computed 
analytically (Lerch, 1986). Then its differential 
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Au = AecosM - eAMsinM 


(4.48) 


is also a non singular variable and can be computed analytically on 
the reference orbit. Considering that a and Aa, are also non 
singular, the first order radial orbit error of (4.46) can be 
approximated by 


Ar G1 = Aa x ' - ‘(aoAei+eoAaJcosMoCt) + a 0 e 0 AM!sinM 0 (t) (4.49) 


which is a non singular expression for the first order radial orbit 
error although the individual errors Ae t and AM, are not. Equation 
(4.49) provides the radial orbit errors of the precise ephemeris, 
computed approximately on the reference orbit. Tests on the 
accuracy of (4.49) are to be made in Section 4.7. 

To compute Aa t , Ae!, AM t the rates Aa, Ae, AM need to be 
integrated, in the sense 


(•At . 

As x (t) = As(t)dt (4.50) 

J o 

where At = t - t„ 

The rates Aa, Ae, AM are computed on the reference orbit using 
equations (4.5) and (4.15), with (4.15) being applied for the potential 
coefficient corrections ACi ma . Dropping the subscript "o" of the mean 
elements for simplicity, we obtain the following expressions for the 
rates 


A&(t) = -l AC. a. sin(^. t + f . ) 

impqa * ma *"»pq *mpqa 0 


(4.51) 


A6(t) = -Z AC. e. sin(V'. t + f. ) 
Impqa £ma £m P t f * m Pq im Pq a o 


(4.52) 


AM(t) = X AC. M. cos(i/ t t + Y'. ) 

£mpqa *«" a £ mpq *«"pq *<"Pq a o 


(4.53) 


where 
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; = -2tf 

^mpq na 2 


« ft) F ^p< i)G 4„ (e > ( ‘- 2l>+<!) 


(4.54) 


e. = fe) F. (i)G. (e)[(l-e»)*(i-2p+q) - (J-2p) ] (1-e 2 )* 

Jtmpq na 3 e k *mp *pq 


(4.55) 


M 


£mpq na 3 


*ft]\.p (1 > [- + 2 < 1+1 > G «p, <e) ] (4 - 56) 


V'. = (i-2p+q)M + (4-2p)u + m(O-0) 


cmpq 


(4.57) 


= (*“2p+q)M 0 + (i-2p)w 0 + m(Q o -0 o ) 


smpqa 0 


- | [a 4 (l-(-l) 1 -")] 


(4.58) 


In the above equations V'lmpq describes the frequency of each 
particular summation term, and Vimpqao gives the corresponding 
phase. Each frequency is a linear combination of the following 
frequency components: the orbital frequency w + M (one cycle per 

revolution), the apsidal frequency 6 of a complete revolution of the 
perigee and the frequency 6-9 of a complete rotation of the earth 
with respect to the precessing orbital plane. One such rotation is 
called a nodal day. 

Equations (4.51) - (4.53) are sinusoidal oscillations with time 

independent amplitudes frequencies and phases. These equations are 
generally valid for a short time period because of the orbital decay 
due to surface forces that can change the overall characteristics of 
the orbit. Of extreme importance is the correct determination of the 
reference orbit since errors in this orbit can substantially distort the 
computed frequencies and phases. On the contrary, the amplitudes 
are not very sensitive to such errors. From equations (4.54) - (4.57) 
and the definition of the eccentricity functions we can conclude that 
the major contribution to the semimajor axis error arises with q = 0 
while for the eccentricity and mean anomaly such contribution is 
given from terms involving q * 0. Since altimetric orbits are almost 
circular, truncation of the series at q = *1 is adequate. 

When the frequencies pq are different from zero each term of 
the summations in (4.51) - (4.53) can be normally integrated as a 
cosine or sine function. The total integrated effect is then 
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Aai (t) = AC, :* mpq cos(f, At + y ) - Aa 0 (4.59) 

" — # mpq *mpqa Q 


A L, * 

«mpqa *ma -a 

**mpq 


Ae^t) = ,I AC. _lm£a C0S (Y- At + V* ) - Ae c 

tmpqa *ma *mpq «mpqa Q 1 

**mpq 


(4.60) 




i£a 


AMl(t) =lLa- AC i- + ) ~ AM 0 (4.61) 


mpqa *ma ^ «mpq *mpqa Q 

^fmpq 


where the summations are over the indices that do not make mpq = 
0. Aa 0 , Ae 0 , AM 0 are the quantities computed at At = 0 to make 
Aai(t 0 ) = Ae & (t 0 ) = AM t (t 0 ) = 0, This is physically meaningful because 
errors in the forcing function cannot affect the orbit instantly. They 
only accumulate with the evolution of time. 

Substitution of (4.59) - (4.61) into (4.49) gives the periodic part 
of the first order approximation of the radial orbit error. After 
applying trigonometric identities to the resulting products of cosines 
and sines we obtain 


Ar{jt(t) =*L.„ A C, ^ A° cosfV'* A t + f 


fmpqa ^ma L ^mpq 


&mp q 


mpqa D 


+ a\ cos(V'. , , .At + jp ] 

*mpq *mp( q+O *mp(q+l)a 0 


-l 

+ A. cos 

*mpq 




where 


- Aa 0 + (e 0 Aa 0 +a n Ae 0 )cosM 0 (t) - a o e o AM o sinM 0 (t) (4.62) 
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wit.h a similar expression for ^£ mp ( q ± l ) a0 . From (4.62) it be >c:oines 
obvious that each component of the disturbing function produces, in 
addition to an error component of frequency ip, an excitation 

modulated by *M (or by 1 cy/rev). The contribution at frequency 
is given by the semimajor axis errors while the modulated components 
are a result of errors in the eccentricity and mean anomaly. So any 
long period effects in these two elements will produce an error with a 
frequency close to 1 cy/rev. On the contrary, errors close to 1 
cy/rev will produce long wavelength errors as well as errors with a 
frequency close to 2 cy/rev. In addition to these periodic terms the 
radial orbit error contains a constant bias and a 1 cy/rev terms. 
These two error components are a function of all the potential 

coefficient terms. This is in contrast to the periodic errors that are 
functions of particular degrees and orders giving rise to 
corresponding frequencies. 

4.5.4 Resonant Effects 

When the disturbing function has a frequency that is identical 
or very close to one of the eigenfrequencies of the homogeneous 

partial differential equations that correspond to the equations of 

motion (4.4), then resonance occurs. Simply stated, resonance occurs 
whenever the satellite is subjected repeatedly to the same force field 
resulting into large very long wavelength oscillations or even a 
secular change in the satellite motion. A secular change is called 

perfect resonance while the long period oscillations are called deep or 
shallow resonances. 

A resonance appears in the form of zero divisors (perfect 

resonance) or small divisors (deep or shallow resonance) in the 

solution of the linear equations of the type (4.51) - (4.53). The 
existence of small divisors though does not imply that the 

perturbations have any infinite or wild behavior, but the fact that 

the analytical solutions break down. In order to see what particular 
conditions create resonant effects, equation (4.57) has to be written 
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as 


■i/i. = (i-2p+q) (m+M) - q« + m(6-0) (4.65) 

* mpq 

This equation becomes zero for all even zonal terms for which q = 
0 and p = t/2. All those terms produce secular changes (perfect 
resonance) in the Keplerian elements and their errors. Long period 
perturbations on the other hand (deep resonance) are produced by 
the odd zonal terms for which p = (f+l)/2 * 1 and q = *1. Then 
fompq = *«• Since w * 1-2 degrees per day these perturbations have 
periods of about 200 days and large amplitudes. Another case in 
which resonance occurs is when the frequencies involved in (4.65) 
become locked together. In other words when the orbital frequency 
(w+M) is related to the nodal frequency (Q-0) by a simple integer 
ratio. When this happens the satellite passes through the same 
geographic region and is subjected to the same gravitational features, 
resulting into a gradual accumulation of the corresponding 
perturbations that have a frequency equal to q«. When either q or w 
are equal to zero then there is a perfect resonance (w = 0 is the 
resonant constraint for frozen altimeter orbits). When qw is different 
from zero then there is a deep resonance. 

In reality « can never be identical to zero, or even be kept very 
close to zero throughout the satellite lifetime because the orbit drifts 
with time. This constraint can only be maintained if small corrections 
are applied to the orbit by firing rockets onboard the spacecraft. 
Furthermore, the orbital and nodal frequencies can never be perfectly 
locked together because of interactions with periodic perturbations 
and the surface force effects. So instead of perfect and deep 
resonances we have the so called shallow resonances which occur 
when the orbital frequency and the nodal frequency are very close in 
forming an integer ratio. As an example of shallow resonances 
consider a satellite with an orbital frequency of 14.3 cy/day. Then 
the first shallow resonance occurs at m = 14 (primary shallow 

resonance for which ^ i mp q » 0.3 cy/day) the second at m = 28-29 
(secondary shallow resonance) and so on. 

Next in importance to perfect, deep and shallow resonances are 
the effects of the terms with J8-2p+q = 0 and m * 0. Then the 
frequency reduces to -q« + m(O-0) which can never become zero since 
all perturbations of importance have a q < 2 in nearly circular orbits. 
This frequency becomes small for small values of m. And since m 
multiplies the nodal frequency the corresponding terms are called 
m-dailies. 

The resonant effects in the Keplerian elements can be seen by 
examining equations (4.51) - (4.56). It is obvious that Aa t is not 
affected by secular and long period effects of zonal origin and by the 
m-dailies because in all cases f-2p+q = 0 and so (4.54) becomes zero. 
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Secular changes occur only in AM X while long period effects of zonal 
origin and m-dailies occur in both Ae x , AM X . Shallow resonances 
occur in all Aa x , Ae x , AM X . Long period effects due to shallow 

resonances and m-dailies can generally be computed in a normal way 
using the theory as described in the previous Section. The secular 
changes in AM X can be computed by direct integration of (4.53) 
considering ^i mpq = 0. Then 

AM r ( t ) = J AC^M, .At (4.66) 


where t is even. For long period perturbations of zonal origin and 
for a small integration period (i.e. about 7-10 days) the oscillations in 
Ae and AM expressed by (4.52) and (4.53) can be approximated by a 
linear trend. Then these equations can be integrated to 


Ae r (t) = ] 

i 

C AC. e. £± i Atcos« 0 

i ioo ° 

2 

(4.67) 

AM r (t) = ] 

j 

\ AC i 0o M i 0 i*t ?1 Atsin "o 
2 

(4.68) 


where i is odd. Note that for a longer integration period a quadratic 
term should be included to account for the non linear effect of the 
long period oscillation. The general form of the equations (4.66), 
(4.67) and (4.68) is 


Ae r (t) = e r At 

(4.69) 

AM r ( t ) = M r At 

(4.70) 


where 


e r = ~ X AC. e. sinf. 

Impqa £ma *<"pq £ n>pqa 0 


(4.71) 


M r = X AC. M. cosV. 

impqa * ma *mpq *mpqa 0 


(4.72) 
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These equations are more general than the previous ones, because 
they can be used to compute any type of resonance that may occur 
and not only the ones of zonal origin. This can be done by 
establishing a threshold while computing fompq and then using the 
above equations for all combinations of A, m, p, q that give rise to 
i^mpq smaller than the threshold. Since a resonance can also occur in 
the semimajor axis a similar equation for Aa r has to be included. Aa r 
as well as Acj r (which is to be used in the next section) are given by 


Aa r (t) = a r At (4.73) 

Aw r (t) = « r At (4.74) 


where 


= - X 


AC, 


cmpqa 


tma mp q 
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*ma *mpq 
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Cmpqa 0 
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Kmpq 



aF lmp( i ) Q 
3i Apo 

tani(l-e 2 )* 


(e) 


(4.77) 


Equation (4.77) can be obtained by substituting the forcing function 
into the equation of motion (4.5) for «. 

Using (4.69), (4.70) and (4.73) in (4.49) the first order radial orbit 
error due to resonant terms is 


Ar G i(t) = a r At - (a 0 e r +e 0 a r ) AtcosM 0 (t) 


+ a 0 e 0 M r AtsinM 0 (t) 


(4.78) 
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which contains a linear trend and a 1 cy/rev component with an 
amplitude linearly increasing with time. This type of error is also 
zero at t = t Q . 


4.5.5 Second Order Radial Orbit Error 

In the previous sections the first order radial orbit error has 
been computed. In these computations it has been assumed that the 
Keplerian elements that are used as arguments in equations (4.59) - 
(4.61) are the mean elements which are time independent and contain 
no errors. In reality these elements are the osculating elements and 
have perturbations that interact with the computed errors to produce 
second order effects. These second order errors in the radial 
distance are 


Ar G2 


3r I - 3§ 

di J 9s 


Ai t dt 


(4.79) 


Since (4.79) involves the interaction of the singular elements e, 
M with their first order errors Ae lt AM 2 the non singular 
transformation (4.47) has to be applied. Then (4.79) becomes 


Ar G2 


3r f as'.- .. 

al'J Is As ‘ dt 


where 


s' = (a,u) = (a.ecosM) 


(4.80) 


(4.81) 


Since the errors A§! are generally small, all the interactions with 
the perturbed elements s' will be negligible with the exception of the 
resonant errors which accumulate with time, and the constant error in 
the semimajor axis. The latter one can be quite important because of 
its interaction with the mean motion n through Kepler’s third law. 
Then, using (4.81) we can write (4.80) as follows 


Ar G2 = -(1-u) } |J Aa 0 dt + a j |jj Aa„dt 


+ (1-u) } |J e r Atdt - a { j* e r Atdt 
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+ (1 -“> I 1 M r Atdt -.J|g S’-Atdt 


+ (1-u) } || « r Atdt - a | « r Atdt (4.82) 

Equation (4.82) is computed in Appendix C. The computation has 
been made separately for the coupling effects of the semimajor axis 
and the variable, u, so that their individual second order effects can 
be determined. These individual effects are given by equations (C.17) 
and (C.37), while their combined effect on the radial distance is given 
by (C.38). All the integrations have been made on the reference 
orbit. 

From (C.17) we observe that the largest second order contribution 
in the semimajor axis is an oscillation with a frequency of 2 cy/rev 
and an amplitude linearly increasing with time. A similar error with a 
smaller amplitude exists in the variable, u. These two errors have a 
different sign so the resulting error in the radial distance is largely 
reduced. Another oscillation with a frequency of 1 cy/rev and an 
amplitude linearly increasing with time is the result of the interaction 
of the constant bias with the mean motion. These two components of 
the radial orbit error are the following 


Ar G2 “ a o e o(| “ Aa 0 ] AtsinM 0 (t) 


- | a 0 J 2 [|*] [| | Aa 0 + M r + w r ] Atsin 2 isin2[M 0 (t) + u 0 (t)] 

(4.83) 


The magnitude of the oscillations in (4.83) is clearly dependent 
on the satellite specifications and on the errors of the gravity^ field 
that is used for the orbit determination. The quantities Aa 0 , M r , Z r 
can become very small when an accurate gravity field is used. 
Additionally since they do not have any periodic behavior they can be 
removed rather easily during the orbital adjustment by comparing the 
computed orbit with the observations. Care should be exercised in 
this type of removal because if the observations are not globally 
distributed and accurate, the removal of these quantities will not be 
efficient. In the context of this study we will assume that no such 
removal is made. 
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4.6 The Radial Orbit Error Due to Initial State Vector Error 

The initial state vector that is computed from the least squares 
fit of the various force models to the satellite observations, contains 
errors arising from errors in the force models (and primarily the 
gravity field), and errors in the station coordinates and tracking to 
the satellite. When the initial state vector is used for the orbit 
integration, its errors propagate into the computed orbit. This error 
propagation in terms of the radial distance is given by equation (4.22) 


lr i = ai 117 15 ' (4 ' 84> 

The analytical determination of Arj is possible along the lines 
developed so far in this Chapter. Using equation (4.28) the 
osculating elements s can be written 


i = s 0 + ij + i 2 


(4.85) 


where s 0 is the vector of the mean elements, s t is given by (4.29) 
and s 2 represents the second order coupling effects between the 
different harmonic terms and is equal to 


s 2 


1 If § ‘ dt 


(4.86) 


Since s t and i 2 are computed using the mean elements s 0 which are 
determined from the initial elements ij, all three terms of (4.85) are 
affected by errors bsj. Then 


Ari 


3r [ 3s„ [ 3s, 1 
3s L 3s j SSj 




(4.87) 


The last partial representing effects of the order of J \ is about three 
orders of magnitude smaller than the other two partials and since As x 
is generally small, this last term is not expected to have any 
significance and therefore it can be dropped. As it is described in 
section (4.5.2) the mean elements are usually computed by removing 
the J 2 perturbations from the initial elements, or by averaging the 
osculating elements, whenever available. In both cases the errors of 
the initial elements propagate directly into the computed mean 
elements. So the partial SSfj/ds! is 
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= I + 0(J 2 ) 


(4.88) 


where I is a unit matrix. Computing also the derivatives 3r/3s on the 
reference orbit we obtain for the first part of (4.87) 


Ar xi = Aa x - (a 0 Ae I +e 0 Aa I )cosM 0 (t) + a 0 e 0 AM x sinM 0 (t) (4.89) 


The errors Ar xi have an identical form to the errors in (4.62) that 
are functions of the constants of integration Aa 0 , Ae 0 and AM 0 . They 
both have a constant part and a 1 cy/rev part but with different 
signs. So initial state vector errors can either reduce or enhance the 
constant and 1 cy/rev errors of gravitational origin. 

The second part of (4.87) can be written as 


* r » ■ H \ Si i 5 i dt < 4 - 9o > 

which is identical to equation (4.79) describing the second order 
effects of gravitational origin. So equation (4.90) can be determined 
in a similar way by introducing the transformation (4.81). Then the 
only initial error that can produce any effect of significance is the 
one in the semimajor axis. Using the equations in Appendix C for Aaj 
we obtain 


Ar I2 = -a 0 e 0 [| ~ Aa x ] AtsinM 0 (t) 

+ | a 0 J 2 (^j) (| Aa x ] Atsin 2 isin2(M 0 (t) + «„(t)) (4.91) 


which again is of the same form with the corresponding error of 
gravitational origin. The interaction of Aa x with the other terms 
other than J 2 are negligible as long as Aa x is small. 


4.7 Validation of the Theory for Errors of Gravitational Origin 

Throughout this Chapter the complete theory for the modeling of 
the radial orbit error has been developed, and formulas providing the 
first order and second order effects of the error have been 
computed. It remains to be seen how accurate this formulation is, so 
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that it can be used for real applications. 

The validation of the theory aims primarily at the assessment of 
the accuracy of the computed first and second order radial error. As 
a by product one can obtain estimates on the accuracy of the 
determination of the reference orbit and on the accuracy in the 
computation of the semimajor axis error Aa and the error aAu of the 
scaled non singular variable au. Furthermore, it can be shown that 
individual errors Ae and AM are indeed non estimable from the theory. 

To make the validation, two numerically integrated orbits were 
provided to us by NASA, Goddard Space Flight Center (Klosko, private 
communication, 1986). These orbits were computed for a six day 
Seasat arc at 2 minute intervals using the GEM10B (Lerch et al., 1981) 
and GEM9 gravity fields and the same initial state vector. Both fields 
were used up to harmonic degree 10 so that to facilitate the 
computations. The satellite state vectors for the two orbits were 
initially in terms of rectangular coordinates and velocities. So they 
have been transformed into Keplerian elements using the 
transformation described in Appendix A. The initial elements used for 
the generation of both orbits are given in Table 1. 


Table 1 

Initial State Vector Used for the GEM10B and GEM9 Numerically 

Integrated Orbits. 


Epoch 

780923 

Semimajor Axis 

7177305.511 m 

Eccentricity 

0.00086 

Inclination 

108T0077 

R.A. of Asc. Node 

160T 9817 

Arg.of Perigee 

o:o 

Mean Anomaly 

o:o 


For the determination of the mean elements the time averages of 
a, e, i as well as the average rate of 0 have been computed using the 
GEM9 ephemeris and are shown in Table 2. Individual mean rates of « 
and M could not be computed because of the erratic behavior of these 
elements which are singular for almost circular orbits. Instead, the 
rate of the total argument « + M, which is a well defined quantity was 
determined and is also shown in Table 2. 

The mean elements have also been computed using the analytic 
formulation. The mean semimajor axis was computed by (4.37) and the 
eccentricity by (4.41). The mean inclination is not computed because 
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of the small variations observed in the ephemeris of the GEM9 orbit. 
For the computation of the rates «, Q, and M, equations (4.35), 
providing the J a secular perturbations, have been used. The 
computation of a has been repeated including the secular and long 
period perturbations from all zonal coefficients with no change in the 
result. The mean elements a, e and rates a, M, « + M, Q are also 
shown in Table 2. 


Table 2 

Mean Elements and Rates Computed by Numerical Averaging and 

Analytical Methods. 


Element 

Averaging 

Analytical 

a 

7168889.259 m 

7168980.740 m 

e 

0.0008 

0.00084 

i 

108.017 

— 

w 

— 

-l!727/day 

M 

— 

5146:699/day 

<u + M 

5144:966/day 

5144!972/day 

Q 

2T053/day 

2!044/day 


From Table 2 we can see the excellent agreement between the two 
methods in all the estimates. The agreement in the mean rate u + M 
indicates that the individual rates are also correct. From Tables 1 
and 2 it turns out that, with the exception of the semimajor axis the 
mean elements and the mean rates can be very well approximated by 
their initial values and the rates computed from J 2 , respectively. 
Note that for the computation of M the mean semimajor axis has to be 
used. 

The GEM10B and GEM9 numerically integrated orbits were 
differenced to generate "errors" in the semimajor axis, the scaled 
variable au, and the radial distance. These errors are shown in 
Figures 1, 3, 5 while their amplitude spectra, computed by a Fast 
Fourier Transform technique, are shown in Figures 2, 4, 6. From 
these figures we can see that the radial orbit error can reach 
magnitudes of about 3-4 meters with the major error being at 1 
cy/rev. Additional energy is found in the frequencies around the 1 
cy/rev frequency primarily due to the m-dailies, while practically no 
energy exists above 3 cy/rev. The 1 cy/rev error comes from aAu 

since there is no 1 cy/rev error in Aa. The RMS error Ar is 1.74 

meters. The error Aa on the other hand has a constant bias of -1.007 
meters and an RMS value of 1.33 meters. As it can be seen from 

Figures 1 and 2, Aa is increasing in time with an apparent 2 cy/rev 

frequency of large amplitude. Additional frequencies as high as 5 
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cy/rev with smaller amplitudes do exist. All these errors are largely 
reduced when Aa is merged with aAu to create the radial orbit error 
Ar. 


The errors Aa, aAu and Ar have also been computed analytically 
using the differences GEM10B - GEM9 up to degree 10 as potential 
coefficient errors. The first order equations (4.59), (4.62) and (4.78) 
have initially been used. First order results indicate that the 
summation part of (4.62) does not give rise to any errors at 
frequencies of 0, 1, and 2 cy/rev. The constant bias and 1 cy/rev 
errors come exclusively from the last terms of (4.62). Additional 1 
cy/rev errors with a time dependent amplitude arise from the long 
period errors in eccentricity which build up to about 70 cm after six 
days. On the contrary the secular and long period effects of mean 
anomaly errors are negligible (3 mm after six days). 

The first remarkable agreement of the first order results with the 
numerical ones is that the constant bias has a value of -1.009 meters 
which differs from its true value by only 2 mm. In a more detailed 
comparison it is seen that the first order theory is able to model all 
the frequencies including the 1 cy/rev terms of constant amplitude 
with a remarkable accuracy. Indeed the discrepancies from the 
numerical results have systematic patterns of 1 and 2 cy/rev with 
time dependent amplitudes, which are characteristic of the second 
order effects as seen from equation (4.83). With these effects 
included, the agreement in the semimajor axis is excellent. The 
maximum discrepancy is 5 cm and the RMS discrepancy is 1.5 cm. 
These discrepancies are shown in Figure 7 while their amplitude 
spectrum is shown in Figure 8. This spectrum shows frequencies at 
1, 2 and 3 cy/rev with a maximum amplitude of 8 mm. 

On the contrary the agreement in terms of the radial orbit error 
is not as good. From Figures 9 and 10 we can see that there is still 
a systematic discrepancy at 1 cy/rev with time dependent amplitude 
and an RMS value of 20 cm. Superimposed is a discrepancy of daily 
frequency and small magnitude. It turns out that the 1 cy/rev 
discrepancy is equal to (n/a)Aa 0 AtsinM which implies that the first 
term of (4.83) should be three times smaller than it is given. This 
systematic discrepancy is of unknown nature and more investigation 
is needed to resolve it. If we account for this correction, the overall 
discrepancy in Ar has an RMS value of 1.2 cm with a maximum value 
of 5 cm. 

The comparison between the two sets of errors has shown that 
the analytic formulation can provide a very accurate modeling of the 
radial orbit error and can be used to form the detailed observation 
equations. Indeed errors in the theory (with the exception of the 
error in (4.83)) are expected to be smaller than 1 cm when the theory 
is applied for future altimetric missions, like the Topex/Poseidon 
mission, since the radial orbit errors are also expected to be quite 
smaller than the simulated errors used in this comparison. 






TIME IN DATS 

Figure 3. Error in the Non Singular Variable au Computed by 
Differencing Two Numerically Integrated Orbits Based on 
GEM9 and GEM10B to Harmonic Degree 10. 
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Figure 4. Amplitude Spectrum of the Error in the Non Singular 
Variable au. 
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Figure 11. Discrepancies in the Computation of the Radial Orbit Error 
After Additional Empirical Correction Has Been Applied. 
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Figure 12. Amplitude Spectrum of the Discrepancies in the 
Computation of the Radial Orbit Error After Additional 
Empirical Correction has been Applied. 


CHAPTER V 


ALTERNATIVE EXPRESSIONS OF THE RADIAL ORBIT ERROR 


5.1 Introduction 

In the previous Chapter the complete formulation leading to 
detailed equations for the orbit error has been developed. These 
equations, given a set of potential coefficient errors and initial state 
vector errors, are able to provide the radial orbit error over a short 
period of time. Closer examination of them reveals some properties of 
the error, like the strong 1 cy/rev term arising from all potential 
coefficient errors, the linear time dependence of the amplitude of the 
error arising from resonant terms and the constant bias and others. 
Computation of individual summation terms can even reveal information 
about the contribution of specific degrees and orders to the error. 
Additionally it can be seen that different combinations of £, m, p, q 
produce oscillations of the same frequency and phase but of different 
amplitudes. Finally the geographic location of the subsatellite point 
at a particular epoch for which the radial error is computed, can be 
found by simply converting the Keplerian elements at epoch to 
spherical coordinates. 

This type of formulation though does not lend itself to a 

systematic analysis of the properties of the orbit error. Furthermore 
its implementation on the computer is quite expensive. So 

transformations of these equations are necessary. One such 
transformation involves the grouping of all terms that produce the 
same frequencies to create a Fourier series expression. Then the 
frequency content of the error can be easily determined. 
Furthermore, the cumulative contributions by degree and order are 
easier to compute. Another formulation that can be obtained involves 
the systematic transformation of the orbit error from the inertial 

coordinate system to an earth fixed system. This formulation allows 

for the determination of the geographic patterns of the error. 

In this Chapter the Fourier series formulation and the geographic 
representation of the orbit error are going to be developed. Then 
some numerical results using a set of "known" potential coefficient 
errors will be given for both methods. Finally some alternative forms 
for the orbit error used by other investigators will be derived. 
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5.2 Fourier Series Formulation 

5.2.1 Determination of the Fourier Frequencies 

The frequency of each individual term contributing to the orbit 
error is defined by equation (4.65) 


= (*-2p+q)(M+£) - q« + m(Q-0) (5.1) 

*mpq 

The integers £, m, p, q can be any combination of integers in the 
ranges defined by 


* = [ 2 ,* max ] 
m = [0,<] 

P = [0,1] 
q = [-1,1] 


(5.2) 


where £ max is the maximum harmonic degree for which the satellite 
shows any sensitivity to the gravity field. The range of q has been 
restricted to only three values for the case of near circular orbits. 

In (5.1) the most . dominant term in characterising the frequency 
content is (£-2p+q) (M+w) or k cy/rev. The second largest term is 
m(fl-0) or about -m/14.3 cy/rev while qd> is very small and can be 
neglected. Considering the range of variation of p and q a certain 
harmonic degree £ will generate frequencies in^ the range [-£-l,£+l] 
cy/rev. Allowing also for the modulation by *M this range becomes 
[-£-2,£+2] giving a total of (2£+5) different frequencies separated by 

1 cy/rev. Around each of these frequencies there is a number of 

other frequencies generated by the terms m(Q-0). So we can view the 
frequency spectrum of the orbit error for a particular degree as 

being composed by (2£+5) distinct frequency bands 1 cy/rev # apart, 
each of them containing m different frequencies separated by (Q-0) or 
1/14.3 cy/rev. The total number of frequencies is (£+l)(2£+5). 

Obviously some frequencies belonging to different bands and to 
different orders can become identical or have a very small separation. 

It turns out that for all harmonic degrees there are always 
combinations of p and q that generate the same frequency. In other 
words harmonic coefficients of different degrees and same order are 
lumped together in the same frequencies. The total number of 
frequencies generated by a gravity field up to £ max is 

(^max+lME^max+S) while the maximum positive and negative frequencies 
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are 


VW = (^max+2) cy/rev 

(5.3) 

= t-(^max + 2) “ * ma x/14.3] cy/rev 

It is of interest and very important to investigate what 
combinations of i, p, q give rise to waves of the same frequency 
band. For m=0 we have three possible combinations of i, p, q that 
provide the same frequency, k. Denoting by i = (-1,0,1) the 

modulations (-M,0,M), these combinations are 


£-2p+q = k 

for i = 0 

(5.4) 

i— 2p+q+l = k 

for i = 1 

(5.5) 

J-2p+q-l = k 

for i = -1 

(5.6) 


All the possible values of p and q for a particular i that satisfy 
these equations are given in Table 3 . 


Table 3 

Combinations Between i, p, q that generate a frequency k. 


modulation 

i-k even 

■8-k odd 

p»q 

P>3 

i = 0 

±± 0 
2 * U 


i = 1 

i+2-k . 

2 * 1 

i-k . 

2 ' 1 

i+i-k 
2 * U 

i = -1 

t- k , 

2 ’ 1 

i-2-k . 

2 

i-i-k 
2 ’ U 
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From Table 3 it can be seen that each degree t has either four 
or five components contributing to the same frequency k, depending 
on whether i-k is odd or even. Obviously the minimum harmonic 
degree A that can give rise to a frequency k is given by 


A 


ml n 


raax( lkl-l,m,2) for i=0 
max( lkl-2,m,2) i=l 

max( Ikl ,m,2) i=-l 


(5.7) 


Using the above Table and equations (5.4) - (5.7) a frequency of km 
cy/rev can be defined and written as 


Vkin = k(M+«) + m(Q-0) 


(5.8) 


5.2.2. Determination of the Fourier Coefficients 

As it is shown in Chapter 4 the radial orbit error is basically 
composed of a constant bias, a 1 cy/rev term with constant and time 
dependent amplitudes, a 2 cy/rev term with time dependent amplitudes 
and a sum of many periodic terms containing all the frequencies as 
defined in Section 5.2.1. The first three terms are trivial to bring 
into a Fourier series form. The summation part of (4.62) can be 
rewritten as 


Ar = X I Ar 

k m 


(5.9) 


where 


Ar 


km 


= z 

i 


n 


z 

p» q 


. cos(V’ At 

&mpq km 


^Ampqal 


(5.10) 


which shows that a wave of frequency km is a sum of waves of the 
same frequency with amplitudes and phases that are generally not 
identical. So in order to form the main wave of frequency km all the 
components with different phases have to be combined. From (4.16) 
we obtain 
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4 = (*-2p+q+i)M 0 + (l-2p)«„ + m(O o -0 o ) + C (5.11) 

^ fnp qa 1 g 

where 


c = “ |(a+|(l-(-l)^-Jn)) 


(5.12) 


Depending on whether 4-k is even or odd, the combinations of 
£,p,q as shown in Table 3 lead to the following different phases that 
are shown in Table 4 for m=0. 


Table 4 

Phases for p,q from Table 3 that generate a frequency k. 


modulation 

i-k even 

<*-k odd 

p,q 

p,q 

i = 0 

k(M 0 +« o ) + C 

k(M 0 +« 0 ) + C - w„ 
k(M 0 +w 0 ) + C - u 0 

i = 1 

k(M 0 +w 0 ) + C - 2«„ 
k(M 0 +w 0 ) + C 

k(M 0 +w 0 ) + C - « 0 

i = -1 

k(M 0 +« o ) + C 
k(M 0 +« 0 ) + C + 2 w 0 

k(M 0 +w 0 ) + C + « Q 


Table 4 shows that the individual phases are all composed by a 
common phase which is shifted by integer multiples of o 0 . Denoting 
by 


f kmo = k(M 0 +« 0 ) + m(Q o -0 o ) 


the phases of equation (5.11) become 


•\h 

*mpqa 1 0 


-if + C + jw 

kmo 


(5.13) 


(5.14) 


where j takes values in the range [-2,2] consistent with Tables 3 and 
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4. After some manipulation equation (5.10) becomes 

4r „„ = ? I 

- B ».< si "« k „ At ^„. +c) ] 

where 

A . = A, . + A . cosw 0 + A . cos2cj„ 

km* km* o km-Gi 0 km^2 ° 

B, * = B, . sin« 0 + B . sin2cj_ 

km* km* i 0 kmJ0 2 ° 

and 


(5.15) 


(5.16) 

(5.17) 


A 


km£ o 


A . 

km^ l 


B , 

kmj 8 1 


^kmj® 2 


^km^ 2 


A «Jt + ‘Ijbi-, + * ‘“A even 
0 t - k odd 

0 

+ A i, „* + |~ k -i + k lJ*±=±o + A f „i * = l= k o 

(5.18) 


0 

~ k °tJ^ -i + h °tJ^- 




O 




“1 
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The coefficients in these equations are obtained from (4.63). In 
equation (5.15) the summation over & starts from i min = 
max( I k I -2,m,2), which is the lower bound of equations (5.7), so that to 
account for the coefficients Obviously the coefficients A°£ mpq 

are zero for * m1n , while the coefficients A -l i mpq are zero for both 
^min (i m i n +l). 

For the final transformation of equation (5.15), to bring it into a 
Fourier series form, we have to use (5.12) for the constant C, write 
the summation over 'a' explicitly, and apply some trigonometric 
identities. Then equation (5.15) becomes 


Ar = C cost!' At + S sin ib At 

km km km km km 


where 


C = D cosoJ' + E sinf 

km km kmo km kmo 


S = E cost!/ - D sinf 

km km kmo km kmo 


with 




km^ 


-E 


km^ 




m even 


i—m odd 


km 


kmi 


km^ 


-m even 


m odd 


and 


D a — AO * A * As * B A 
km* km^ km>fl 


E kmi = AS im A kmi! + AC im B kml 


(5.19) 


(5.20) 

(5.21) 


(5.22) 


(5.23) 


(5.24) 

(5.25) 


Equation (5.18) can also be written in terms of its amplitude and 
phase 


60 


Ar = A cos (if/ At-Y 1 ) 

km km km km 


(5.26) 


where the amplitude at frequency km is 


A = (C 2 + S 2 ) % 

km km km 


and the phase is 


i> = tan' 1 

km L km 


(5.27) 


(5.28) 


Then the total radial error representing the summation part of (4.62) 
can be written in either of the two equivalent forms 


Ar = I I (C cosY' At+S sin At) (5.29) 

k m km km km km 


or 

Ar = E I A cos At-^ ) (5.30) 

k m km 

The bias and 1 cy/rev terms in (4.62) become 


Ar 0 = ~Aa 0 (5.31) 

Ar! = CjCosMAt + S^inMAt (5.32) 

or Ar ! = AjCOsfMAt-M) 
where 

C x = (e 0 Aa 0 +a 0 Ae 0 )cosM 0 + e 0 a 0 AM 0 sinM 0 (5.23) 

Si - — ( e 0 Au u +a C) Ae t , ) sinM 0 + e 0 a 0 AM 0 cosM 0 (5.34) 
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Ai = (C?+S?) Vl (5.35) 

M = tan -1 ^ (5.36) 

V' 1 


while similar expressions exist for the oscillations with amplitudes that 
increase with time. 

Equations (5.29) or (5.30) are very convenient to compute various 
statistics of the periodic orbit error. More specifically the amplitude 
spectrum is readily obtained by the amplitudes A km while the average 
power by degree, order, or cumulative, over the period of the orbit 
can be easily computed. The average power for a Fourier series in 
general, is defined by 


P = 


1 

T 


[l ( a * cosx j t+b , s inx j t ) ] dt 


(5.37) 


where T is the integration limit and is chosen such that it is longer 
than the period corresponding to any x^Xj that is formed after the 
expansion of the square of the series. Integration of (5.37) leads to 


P = | l (a’+b-) 


(5.38) 


Based on (5.38) the RMS orbit error over a period longer than the 
period corresponding to any frequency j will be 


I I Ar I \ % 



(5.39) 


Restricting the summation to be only over k, we can have the 
RMS orbit error by order 

" ir "« = [2 l < c ^ + C)] S (5.40) 

Finally the RMS orbit error by degree, after some manipulation, 
becomes 
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"AH I* = [| ) (^c| m +ASf m )]^ 


(5.41) 


By breaking (5.41) even further, one can obtain the RMS effect for 
any individual degree and order. This effect is summed only over k 
that has a k m)n = t and takes both positive and negative values 
absolutely greater or equal to k mjn . 


5.2.3. Numerical Results 

In this section, some numerical estimates of the radial error, its 
frequency content and its statistics will be given. These numerical 
estimates will be computed based on the orbit characteristics of the 
Seasat satellite and a set of "known" potential coefficient errors up to 
harmonic degree 36. This harmonic degree is chosen because the 
sensitivity of Seasat to gravity field perturbations of higher 
harmonics is negligible, with the possible exception of some resonant 
terms. 

To create the potential coefficient errors the half differences 
between the GEM10B and GRIM3L1 (Reigber et al., 1985) gravity models 
were used. Both these fields are given complete to degree and order 
36 and have been independently derived, although not from entirely 
independent data. They are both produced by combining satellite 
observations, surface gravity measurements and altimeter data. Their 
greatest differences are, the larger number of satellites used in the 
GEM10B solution and the more recent data used in the GRIM3L1 
solution, particularly the laser tracking data of Lageos and the Seasat 
altimeter data. Both of these fields should represent reasonable 

global models of the earth’s gravity field, so their differences, 
divided by two, could be a possible realization of their errors. 

The radial orbit error has been computed for the first six day 
arc of Seasat satellite, using equations (4.62), (4.78) and (4.83) and is 
shown in Figure 13. It may be seen that the most dominant features 
of the error are, the 1 cy/rev effects, as well as the daily terms. 
Note also the gradual increase of the error with time, that is due 
primarily to the eccentricity resonant errors and second order 
effects. 

The error for the same six day arc has also been computed using 
the Fourier series equations (5.16) - (5.29). The agreement between 
the two methods is at the 1 mm level, indicating that the Fourier 
series is equally accurate to the original error expression. After this 
consistency check, the Fourier series equations were used to generate 
the radial orbit errors for all six day arcs of Seasat. The results 
indicate that the periodic part of the first order error (equation 5.29) 
is practically the same for all the arcs. Its variations from arc to arc 
is on the order of 1-2 cm in individual values. The RMS error 
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o 

o 



TIME IN DAYS 


Figure 13. Seasat Radial Orbit Error Computed Analytically Using the 
GEM10B - GRIM3L1 Half Differences as Potential Coefficient 
Errors. 


for any of the arcs is 1.65 meters as computed by (5.39). The 
amplitude spectrum of the error, up to 6 cy/rev, for the first Seasat 
arc (and for all arcs) is shown in Figure 14. In this Figure we 
observe that most of the energy of the error is below 2 cy/rev. 
Tests have shown that the total energy of the error above 6 cy/rev 
is about 5 mm. It can be seen that there is no energy at the 1 
cy/rev frequency arising from implemetation of (5.29). The same is 
true for the zero frequency and the 2 cy/rev frequency. This is 
partly due to the fact that the zonal terms (that give rise to these 
frequencies) are very small, and partly because the sensitivity of the 
spectrum at that frequency is extremely small. This sensitivity has 
been computed by assuming a uniform potential coefficient error, of 
5*10 -9 , for all degrees and orders. This spectrum is shown in Figure 
15, where indeed there is almost no radial error at those frequencies. 

The constant term and the 1 cy/rev error have been computed by 
(5.31) - (5.36) and are clearly dependent on the particular arc for 
which they are computed. The magnitudes of these two types of 
errors are shown in Table 5 for all the arcs of Seasat. As it can be 
seen, both errors have significant variations from arc to arc. The 1 
cy/rev error is the dominant error for all arcs. Using these values 
and the RMS periodic error of 1.65 m, the total RMS error for each 
arc has been computed and is shown in Table 5. 
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14. Amplitude Spectrum of Seasat Radial Orbit Error Using 
the GEM10B - GRIM3L1 Half Differences as Potential 
Coefficient Errors. 
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Figure 15. Amplitude Spectrum of Seasat Radial Orbit Error Using a 
Uniform Potential Coefficient Error of 5-10 -9 . 
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Table 5 


Constant, 1 cy/rev and Total RMS Radial Orbit Error for Seasat Arcs 
Using the GEM10B - GRIM3L1 Half Differences as Potential Coefficient 

Errors. 


Seasat Arc 

Constant 

1 cy/rev 

Total RMS 

1 

0.33 m 

0.79 n 

1.78 m 

2 

-0.44 

1.58 

2.05 

3 

0.57 

2.37 

2.43 

4 

-0.43 

2.96 

2.70 

5 

-0.18 

3.35 

2.90 

6 

-0.51 

2.05 

2.27 

7 

-0.07 

1.30 

1.90 

8 

-0.16 

1.03 

1.82 

9 

-0.47 

2.59 

2.52 

10 

0.75 

2.68 

2.63 

11 

1.10 

3.63 

3.25 

12 

0.09 

2.68 

2.52 


Using equations (5.40) and (5.41) the RMS values by harmonic 
degree, harmonic order and individual degrees and orders have also 
been computed. The RMS errors by order for the (GEM10B - 
GRIM3L1J/2 errors and the uniform errors are shown in Figures 16 
and 17. From both Figures we can see the higher sensitivity of the 
errors to harmonic order 1 (daily errors), to orders 14, 15, and 16 
(primary shallow resonance) and to order 29 (secondary shallow 
resonance). Comparing the two Figures we can see the definite large 
discrepancies of the GEM10B and GRIM3L1 gravity fields between 
orders 5 to 12. Additional conclusions can be obtained from the RMS 
errors by degree, by comparing Figures 18 and 19. More specifically, 
there is a much higher sensitivity of the error to odd degrees, for 
degrees below the primary shallow resonance, and there is a change 
in sensitivity for degrees between resonant degrees. Again it 
becomes clear that the GEM10B and GRIM3L1 fields have large 
discrepancies in degrees 5 to 12. Table 6 containing the RMS errors 
for individual degrees and orders shows basically the same type of 
information. 




RADIAL ERROR IN 


U .00 6.00 12.00 18.00 24.00 30.00 36.00 

HARMONIC ORDER 

Figure 16. RMS Seasat Radial Orbit Error by Harmonic Order Using 
the GEM10B-GRIM3L1 Half Differences as Potential 
Coefficient Errors. 
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Figure 17. RMS Seasat Radial Orbit Error by Harmonic Order Using a 
Uniform Potential Coefficient Error of 5-10 -9 . 
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Figure 18. RMS Seasat Radial Orbit Error by Harmonic Degree Using 
the GEM10B-GRIM3L1 Half Differences as Potential 
Coefficient Errors. 
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Figure 19. RMS Seasat Radial Orbit Error by Harmonic Degree Using 
a Uniform Potential coefficient Error of 5-10 -9 . 
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Table 6 


RMS of Seasat Radial Orbit Error for Each Coefficient Degree and Order 
Using the GEM10B - GRIM3L1 Half Differences as Potential Coefficient 







Errors. 

Units 
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Centimeters 
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5.3 Approximate Expressions of the Radial Error for a Near Circular 
Orbit 

The equations for the radial orbit error that have been derived 
in Chapter 4 and in the previous sections are general equations 
involving no approximations except for the truncation of the series 
containing the eccentricity functions to q = *1 in the Fourier series 
form. The general expression for the Fourier series for any order of 
eccentricity can be easily computed by extending the Tables 3 and 4 
to account for the higher orders. Then equations (5.16), (5.17), and 

(5.18) can also be extended accordingly. 

In this section simpler forms for equations (4.62) and (5.29) will 
be obtained, using the explicit values of G^p^(e) as given by equation 

(4.18) for small eccentricities, and recognizing that Aa is primarily 
sensitive to terms containing Gf po (e) while Ae, AM are sensitive to 
terms containing G^pt^e). A notation consistent with equation (4.13) 
will be adopted. Based on the above we obtain 


f a„) 2(i-2p) „ 

Aa. = na — 6 F. — ? • LL S. 

*mpo ^ a * *mp f0 *mpo 

**mpo 


(5.42) 


Ae. ± =n(^)V + S fi 1 (5. 

imp*! v a > *mp [ 2i) *mpi 2 jb *mp“lj 

**mpi *mp~i 


43) 


AM. = - ^ 

imp*! e t a j imp 


-i+4p+l g + 3i-4p+l g 


I 2 tf. imp-1 2-lfr. * m P 

r imp — 1 r impi 


imp 1 


(5.44) 


In these equations Sf mpq is 


S. = S 

*mpq 


Smpq 


(At) - 


S. (o) 

tmpq 


(5.45) 


and has the functional form given by (4.14). Si mpq is the integrated 
Si m pq with respect to its argument, it has the form 
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sin[ (i-2p+q)M + (i-2p)« + m(Q-0)] 


AS im 


cos[(i-2p+q)M + (i-2p)« + m(Q-0)] 


and is equal to 


(5.46) 


S. = S. (At) - S. (o) 

*mpq *mpq *mpq 


(5.47) 


Using (5.42) - (5.44) in (4.46) and neglecting the term eAa we 
obtain 


Ar . = anf **] F. 

imp t a J imp 


■ ^mpo 




mpo 


+ 4g . 3f 1 fs. cosM + S. sinMl 

2*0 v Impi *mpi > 

**mp l 

+ 4 ?~^ + l fs. cosM - S. sinMl] (5.48) 

t imp-i imp-i J 

'imp— l 

Then using (5.45) and (5.47) and the equalities 
S. (At) cosM + S. (At)sinM = S. (At) 

impi *mpi *mpo 

S. (At) cosM + S. (At)sinM = S 4 (At) (5.49) 

imp-i imp-i impo 

S. (o)cosM + S. (o)sinM = S. (o)cosMAt + S. (o)sinMAt 

*mpi * mp i impo impo 

S. (o)cosM - S. (o)sinM = S. (o)cosMAt — S. (o)sinMAt 

imp-1 imp-1 impo impo 
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we obtain for the periodic first order radial orbit error the following 
expression 


Ar Gi = Ar p - Ar 0 - Ar x 


(5.50) 


with 


Ar D = Z (A. + B. + C. )S 4 (At) 

P imp ^ ^ m P *mp impJ impo 


Ar 0 = 1 A. S 4 (o) 


imp £m P * m P° 


(5.51) 

(5.52) 


Ar x = I r(B + C. )S # (o)cosMAt + (B. - C. )S. (o)sinMAtl 

p L *mp ^mp *mpo *mp *mp «mpo J 


mp *mp *mpt 


(5.53) 


where 


A =J MizM 

*mp *mp ^ 

r impo 


B = J, - 

*mp ^mp Ojb 

*mp l 


c =J 
-tmp *mp 2 t A 

^imp— i 


j < mp = ™( f‘) lr i. P ( i) 


(5.54) 


Equations (5.50) - (5.54) provide an expression for the radial 
orbit error that is equivalent to (4.62) when the eccentricity is small. 
These equations have also been computed by Rosborough (1986). 

The Fourier series representation of Ar p can be easily obtained. 
First the summation of the constant quantities A,j mp , Bi mp , Ci mp is 
rearranged by setting k = i - 2p and summing over m, k. The sum 
of these quantities for any i, m, p is 
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A. + 

*mp 


^imp ^mp 


J^mp 


2 k 


l+2k+l + l-2k+l ~ 
2('yH-M) 2(V^-M). 


which is equal to 



with p - % 
M 


(5.55) 


(5.56) 


Using (5.54) for J£ mp and (4.14) for Si mpo and separating the 
individual phases from the corresponding frequency arguments in the 
cosines and sines we obtain 


Ar. 


^max 

■ .u 


t 


i « 


max ‘max 


t 2 H. fC. cos y At + S. sintf Atl (5.57) 
m =o t = t . ^kml A km r km ikm r km J 

min 


where 


H 


i km 


= a 


ffJ J F t Milihgk 
la J 'iJtJl 0 ( 0 2 -l) 


£ km 


' AC i 
-AS, 


cosy 


kmo 


AS i 

AC, 


sin'V' 


kmo 


(5.58) 

(5.59) 


S ikm 


AS i 

AC, 


cosy 

r kmo 


AC 

AS 


^m 

im- 


siny 

r kmo 


(5.60) 


The summation in £ is made over the harmonic degrees that have the 
same parity with k. The minimum harmonic degree is given by 


i m1n = max(lkl, 2, m+mod( lkl-m,2) ) 


(5.61) 
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Equation (5.57) is identical to the one given by Wagner (1985). 
This equation becomes singular when i& km =0 due to resonance with 
harmonics of order m*0 or when cy/rev due to secular and long 

period effects in the mean anomaly and eccentricity arising from even 
and odd zonal errors. As it is mentioned in Chapter 4 the first type 
of resonance is not likely to occur for a short arc while the second 
type of resonance is computed from equations (4.66), (4.67), (4.68) and 
(4.78). Expressions for these resonant terms, more in the line of this 
section can be obtained by simply using the specific values of e r and 
M r contained in these equations, and the eccentricity functions as 
given by (4.18). For the particular case of near zero eccentricities 
long period effects in the mean anomaly cannot be computed since 
this formulation breaks down. The resonant effects eAM r (t) though 
can be obtained. Making these substitutions we have 


Ae r (t) 



& 

F (i-l)AC. cos w 0 

*o— -r- *° 


(5.62) 


eAM r (t) = 7 X 

* odd 


-pi ne 
*even 



1-1) AC, 
(l+l)«-4) 


sinw n 

o 



(5.63) 


For the computation of the second term of eAM r (t) the approximate 
equality 



» ami 
2 


e 


(5.64) 


has been used. From equation (5.63) it becomes obvious that the odd 
zonal errors have a much greater effect than the even zonal errors, 
since each of the latter is multiplied by the eccentricity. 


5.4 Geographic Representation of the Radial Orbit Error 

The equations that have been developed so far describe the 
temporal characteristics of the radial orbit error in an inertial 
coordinate system. These results are important for understanding the 
nature of the orbit error and its evolution in time as well as for 
computing appropriate statistics. Of equal importance is the 
understanding of the spatial behavior of the orbit error, both along 
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and across the subsatellite tracks. A systematic analysis of the 
spatial behavior of the error requires the functional transformation of 
the Keplerian elements to geographic coordinates in an earth fixed 
system. Such transformation can be viewed as the inverse Lagrange 
transformation where now the S£ mpq functions have to be transformed 
into functions of geographic dependence. In order to do that the 
orbital geometry of the satellite has to be considered. From Kaula 
(1966, p.32) we have the following relationships between Keplerian 
elements and the geocentric coordinates of the subsatellite points. 


0-0 = \ 
cos(a-O) 


(a-0) 

COS(w+f ) 

cos<fr 


(5.65) 

(5.66) 


sin(a-fl) 


sin(w+f )cosi 
cos4> 


(5.67) 


sin(w+f) 


sin$ 

sini 


(5.68) 


COs(w+f ) 



• 2 ^^ 
sin z <M 

sin 2 iJ 


(5.69) 


where a is the right ascension of the subsatellite point. In equation 
(5.69), the plus sign corresponds to an ascending arc (w+f is in the 
first or fourth quadrant) and the minus sign corresponds to a 
descending arc (w+f is in the second or third quadrant). 


5.4.1. First Order Periodic Radial Orbit Error 

For the implementation of the above equations in one of the 
expressions developed so far for the first order periodic radial orbit 
error, the true anomaly f has to be approximated by the mean 
anomaly M. This is acceptable for nearly circular orbits like the 
altimetric ones, where the eccentricity can be assumed to be zero. 
Consistent with this approximation, the expression of the first order 
periodic radial orbit error that has to be used is the approximation 
that is given by equations (5.50) - (5.54). As it will be seen the 
desired transformation of these equations becomes exact. The same 
cannot be true when higher orders of eccentricity are also considered 
because of the additional term -qw that enters the S^ mpq functions. 
Then the radial orbit error of higher orders becomes time dependent 
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without a strict relation to geographic location. 

The quantity that has to be transformed in equation (5.51) is 
SjCmpo* Using (5.65) we obtain 


smpo 


AC i 

-AS, 


m 

m 


cos[ (i-2p) (M+cj) + mX - m(a-Q)] 


+ 


AS i 

AC, 


sin[ (^-2p) (M+u) + mX - m(a-O)] 


(5.70) 


Applying trigonometric identities and setting for simplicity 


£-2p = k 


M+w = v 


(5.71) 


a-Q = u 


we obtain 


Bmpo 


rf ac, ] 


[AS ] 


*m 

-as 

cos(kV"inu) + 

*m 
AC » 

sin(kv-mu) 


cosmX 


+ 


AS i 
AC . 


cos(kv-mu) 



m 

m. 


sin(kv-mu) 


sinmX 


(5.72) 


The next step in the transformation is to express the cosine and 
sine of multiple arcs with respect to single arcs. This computation is 
made in Appendix D. Then 


_ 'iS 1 5? Mkl] (ml. lkl-s . s m-t . t 

cos(kv-mu) - i E I . (-1) cos v sin v cos u sin u 

s— 0 t=0' s 1 

(5.73) 


where s,t have the same parity and 
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3s+t 



2 


k>0 

k<0 


sin(kv-mu) = X I f 1 ^ 1 ) (“)(-l)°cos 

S=0 t=O l S ) It.* 


Ikl 


■S . s 

v sin v 


where s,t have different parities and 


3s+t+l 

2 

s+t+1 

2 


k>0 

k<0 


Using (5.66) - (5.69) in (5.73) and (5.75) we obtain 
cos(kv-mu) = I I ('*') (") (-»" 

I L, | + + 

s+t 


! k 1 +m 

1 y y 

f'ki] 


cos m 4> % Z 

1 s i 

sin 2 <M 

sin 2 iJ 

1 k 1 +m — s — t 
2 

sin<M 

siniJ 


sin(kv-mu) = (*1) 


lkl+m+i 1 


cos m $ 

I k I +m _ s— t 


| l ('s') (t) (- 1 )' 


sin 2 4>l 

2 

fsin+i 

sin 2 iJ 


Uinii 


s+t 


where 0, c are given by (5.74) and (5.76) and s,t 
parity in (5.77) and a different parity in (5.78). In 
form we write 


I kl+m 


Q, (*) = 

*mp 


(±1)* m Q, (♦) 

*mp 


(5.74) 

m-t . t 

cos u sin u 

(5.75) 

(5.76) 


(5.77) 


(5.78) 

have the same 
a more compact 


cos(kv-mu) = (*1) 


(5.79) 
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sin(kv-mu) = (±l) lkl+m R (4>) = (±1)* m+1 R (4>) (5.80) 

*mp *mp 


where the exponents of (±1) have been changed based on the fact 
that k and l have always the same parity. Obviously Qi mp is 
symmetric and Rjg mp is antisymmetric with respect to the equator. 
Substituting these quantities in (5.72) and taking the appropriate 
potential coefficient errors, depending on whether i-m is even or odd, 
we get 


Kmpo 


= Q. (AC. cosmX 

*mp 


AS. sinmX) 

* m 


~(*R. )(AC. sinmX - AS. cosmX) 

* in p » m * m 

when .l-m is even, and 


(5.81) 


S. = R* 

*mpo *mp 


* Q/ 


&mp 


when is odd. 


(AC. cosmX + AS. sinmX) 

*m *m 

(AC. sinmX - AS, cosmX) 

Mm * m 

Combining (5.81) and (5.82) we finally obtain 


(5.82) 


&mpo 


= Q. (AC. cosmX + AS. sinmX) 

*mp Mm Mm 

* R. (AC. sinmX - AS. cosmX) 

*mp Mm Mm 


where 



^■Amp 

R. 

*mp 




even 


(5.83) 


odd 


(5.84) 
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R. 

*mp 

^<£mp 






even 


odd 


(5,85) 


Then Ar p of equation (5.51) becomes 

Ar p (*,X) = I [A. + B, + C } [Q, (AC. cosmA + AS. sinmXl 

r imp ' ^ m P * °ip *mp> L jempv *m ) 

* ~ AS^COSmA)] (5.86) 


Consistent with equation (5.69), the plus sign corresponds to an 
ascending arc and the minus sign to a descending arc. So it becomes 
obvious that Ar p is composed of a mean error that is common to both 
ascending and descending arcs and by a variable error of the same 
magnitude and different sign depending on the direction of the arc. 
Furthermore there is no dependence of Ar p on the actual epoch of the 
orbit integration and on the initial position of the satellite. Both 
components of Ar p are only dependent on the geocentric ctpordinates 
of the subsatellite tracks and the elements a of i„, w, M and 0 that are 
used to compute the quantities A* mp , Bf mp , C* mp , Q* mp and R* mp . 
Since these elements are not expected to change appreciably between 
weekly arcs, the error Ar p is fixed in space throughout the satellite 
lifetime. The spatial variation in Ar p is specified by the behavior of 
the functions Qi mp , Ri mp and the cosines and sines of the longitudes. 
The same property has been numerically observed in equations (4.62) 
and (5.29). 

The constant term Ar 0 of equation (5.52) can be written 



+ JR» (+o)( AC * sinmX 0 -AS. cosmA 0 )l 

*mp *m *m J 


(5.87) 


where (<fr 0 ,A 0 ) are the geocentric coordinates of the subsatellite point 
at the beginning of the arc, and j takes the values 1 and -1 
depending whether the arc starts ascending or descending. 

To transform Ar t of equation (5.53), we need to use the equalities 
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S. (o)cosMAt + S. (o)sinMAt = S,.^ v (o)cos(M+«) 

«mpo "mpo Jmpo 


+ S (A + i) n ,po ( ° )sin(M+Q) 


S. (o)cosMAt - S. (o)sinMAt =8.. . (o)cos(M+«) 

*mpo *mpo (*“ 1 jmpo 


(5.88) 


- S,. (o)sin(M+«) 

1 jmpo 

where the approximation cjAt=0 has been used, to obtain 

cosMAt s cos(cj+M-« 0 -M 0 ) (5,89) 

and similarly for sinMAt* S i mpo can be transformed similar to S^ mpo 
to give 

S* = Q* (AC* sinmX “AS* cosmX) 

*mpo *mp Zm *m 

-( ± R. )(AC. cosmX + AS. sinmX) (5.90) 

*mp Zm Z m 

Then, from (5.53), (5.68), (5.69), (5.83) and (5.90) we obtain 
Art(*,<t» 0 ,X 0 ) = 1 fW* (4> 0 )(AC. cosmX 0 + AS. sinmX Q ) 

S 1X11 £ m p L Imp Zm u Z m 

+ jWj (* 0 )(AC. sinmX 0 - AS. cosmX 0 )l 

*mp zm Zm J 

* f 1- elnS ' f ) [ w « (4*0 )( AC » sinmX „ - AS. cosmX 0 ) 

t Sin lJ imp L * m P *m 

+ jw! (4> 0 )(AC. cosmX 0 + AS. sinmX 0 )l 

*mp zm J 


(5.91) 
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with 


W imp (4,o) " B tmp R (i+i)mp (l,>o) + C imp R (*-i )rap ( * o) 


W,_ _(♦„) = B,__ (♦„) - C. Q,,_ .x (♦„) 


^ mp 


imp (i+l)mp 0 imp (i — l )mp 


W *- ( * o) B Amp R (A+i)mp ( * o) + C imp R ( A-l ) mp ( * o) 


imp 


W Amp (4 * o) B imp Q (A+i)mp ( * o) + C J mp Q (A-l )mp (<,,o) 


(5.92) 


Equations (5.87) and (5.91) give the geographic representation of 
the constant bias and the 1 cy/rev error. Ar„ is constant spatially 
while Ar x is independent of longitude but has a variation with 
latitude. It is composed of a mean error and a variable error. The 
mean error becomes maximum at the nothernmost and southernmost 
latitudes where the variable error becomes zero. On the contrary the 
mean error becomes zero at the equator, where the variable error is 
maximum. 

Both the constant bias and the 1 cy/rev error are dependent on 
the geocentric coordinates of the subsatellite point at the beginning 
of the arc, computed from M 0 , 0 o , and 0 O . These quantities are 

expected to vary appreciably between arcs, so for every arc these 
two types of errors have different magnitudes. Superimposed on the 
arc independent error Ar p they provide the total first order periodic 
error for each arc. The total error is very close to zero in the 
region surrounding the initial point (<t> 0 , A 0 ), and the maximum error 
is in the region with longitude difference of about 180*. 

The geographic representation of the first order periodic error 
due to coefficient errors given by the half differences GEM10B - 
GRIM3L1 has been computed using (5.86), (5.87) and (5.91) for all 

Seasat arcs. Implementation of these equations on a global scale is 
quite expensive in computer time, primarily because the latitude 
dependent functions Qf mp , Rf mp and the equations themselves have to 
be computed for every subsatellite point. An alternative way is to 
compute the error on a regular grid and then to interpolate to any 
subsatellite point to compute the actual error. It turns out that 
computation on a 1* grid can provide interpolated estimates of the 
error that are quite reliable. More specifically the discrepancies 
between the interpolated estimates and the corresponding values 
computed by (4.62) have a maximum value of about 10 cm in the 
northern and southern latitudes, where the error variation is quite 
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rapid, and are generally on the order of 1-2 cm in the equatorial 
regions. 

The error for the first Seasat arc is shown in Figures 20 and 21. 
Figure 20 shows the error in the ascending arcs while Figure 21 
shows the error in the descending arcs. Both Figures have been 
created using a 5* grid with a contour interval of 50 cm. Examining 
these Figures we can observe that the spatial variation of the radial 
orbit error is quite complex. Consistent with the results from the 
analysis of the temporal characteristics of the error, we observe a 
long wavelength along track variation. The across track variation 
though is quite erratic and is by no means of long wavelength 
nature. This is obviously due to the fact that spatially adjacent 
tracks are substantially separated in time and so they have different 
magnitudes in their errors. In some geographic regions this across 
track variability is quite smooth while in others there is a rapid 
change of even 7 meters within distances of 1000 km. The minimum 
error in the ascending arcs is in the equatorial region close to India 
and is very close to zero. This is expected, because the initial 
integration point of the arc is in that region (♦ = -11.97, A = 101.85) 
and so the error should indeed be zero. The maximum error is in the 
region west of the United States with a magnitude of about -5 meters. 
On the contrary, the descending arcs have a maximum error of 4 
meters in the region of the initial integration point and a minimum 
error, very close to zero, in the region of the maximum ascending 
error. Similar patterns of the error, with more or less the same 
magnitude can be observed in maps of the errors in all the other 
arcs of Seasat. The minima and maxima are of course shifted 
depending on where the initial subsatellite point for each arc is 
located, but the overall characteristics are the same. 

5.4.2. Resonant, Second Order and Initial State Vector Induced Error 

To complete the representation of the radial orbit error in terms 
of geographic coordinates the 1 cy/rev and 2 cy/rev components from 
resonant and second order effects as well as the initial state vector 
errors, have to be expressed with respect to (♦,A). 

Starting with (4.89) for Arn we transform cosM and sinM as 
follows 


cosM = cos(M+w) (cosw 0 - «Atsinco 0 ) + sin(M+«) (sinw 0 + «Atcoscj 0 ) 

(5.93) 

sinM = -cos(M+ w ) (sin« 0 + uAtcosw 0 ) + sin(M+w) (coscj q - «Atsin« 0 ) 


(5.94) 
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Figure 20. Radial Orbit Error of Ascending Seasat Arcs Using the 
GEM10B-GRIM3L1 Half Differences as Potential Coefficient 
Errors. C.I. = 50 cm. 
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Figure 21. Radial Orbit Error of Descending Seasat Arcs Using the 
GEM10B-GRIM3L1 Half Differences as Potential Coefficient 
Errors. C.I. = 50 cm. 
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Since a is generally small and the initial state vector errors are 
also small, the time dependent terms can be dropped. Then, using 
(5.68) and (5.69) we have 


Ar 


n 


Aa. 


sin4> 

sini 


(-aAejSin« 0 + aeAMjCosw 0 ) 


* (l- (-aAejCos« 0 - aeAMjSinw 0 ) 


(5.95) 


which demonstrates that the initial state first order error has a 
behavior identical to the constant bias and the 1 cy/rev error of 
gravitational origin. 

For J,he resonant terms, we can use (5.62) and (5.63) in (4.78) and 
neglect a r , to obtain 


Ar r (t) = A t Atcos(« 0 +M) + A 2 AtsinM 

G 1 

where 


A x 


A 2 



(i-l)AC 

0-“T“ 


O 


1 (W)(t-4) 

‘4 2 


AC. 


O 


Using (5.92) and 

cos(« 0 +M) = cos(cj+M) + 
we obtain 

Ar^ (t) = (Ai + A 2 sin« 0 

u 1 


«Atsin(u+M) 


- A 2 cosw 0 d>At) Atcos(w+M) 


(5.96) 


(5.97) 


(5.98) 


(5.99) 



+ (Ai&jAt - A 2 cosw 0 + A 2 sin« Q wAt) Atsin(u+M) 


(5.100) 


where now the terms in wAt 2 are not necessarily negligible, depending 
on the magnitudes of Ai and A 2 . Using (5.68) and (5.69) 


Ar^ (4»,t) 


sin4 

sini 


(AiwAt - A 2 cosw 0 + A 2 sin« Q «At)At 


± 



sin 


sin 




+ A 2 sin« 0 


A 2 cos« 0 «At) At 


(5.101) 


So the resonant induced errors are not strictly a function of 
geocentric coordinates (♦,A). They also depend on the time elapsed 
from the beginning of the arc, in a non linear fashion. They are 
independent of longitude and they contain a mean error and a 
variable error. Unless the coefficients A, and A 2 are very small 
(which imply extremely accurate zonal coefficients) they can reach 
appreciable magnitudes with the evolution of time. Examining Ai and 
A 2 in (5.97), (5.98) it is expected that the major contribution of the 
error will be from the terms containing Ai (resonant errors in 
eccentricity due to odd zonal errors). These errors depend on the 
initial position through 

The second order effects of gravitational and initial state origin 
can be computed in an almost identical way. The 1 cy/rev term of 
Ar G2 in (4.83) as well as the corresponding term of Ar I2 in (4.91) can 
be easily transformed by setting A, of (5.97) equal to zero and A 2 to 
be equal to A 2 where 


A 2 


| {■; (Aai - Aa„) 


(5.102) 


and evaluating (5.101) with A! = 0 and k' 2 


1 ci nd) ^ ' 

Ar ($,t) = — : — r (-A 2 C 0 S« o + A 2 sin« 0 «At) At 
Ga,i2 v ’ ' sini ° ° 


* ^1- ^j^ 2 rj (A 2 sin« 0 - A 2 cos« 0 «At) At (5.103) 
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The final transformation involves the 2 cy/rev terms of (4.83) and 
(4.91). This is a straightforward transformation giving 


Ar 


2 

G2, 12 


(♦,X) 



sin 2 4 )^ sin4> 
sin 2 iJ sini 


A 3 At 


(5.104) 


where 


A 3 


-aJ 2 ( (| -■ ( Aa 0 - Aaj) + M r + « r ]sin 2 i 


(5.105) 


This type of error is also dependent only on latitude and on time in a 
linear fashion. It is independent of the initial position and it has 
only a variable component which becomes zero both at the northern 
and southern latitudes as well as at the equator. 


5.3.3. Observability of the Radial Orbit Error from Crossover 
Discrepancies 

Having been able to understand the spatial behavior of the radial 
orbit error, it is now straightforward to assess its observability 
from crossover discrepancies. This assessment is quite important, 
because, as it is mentioned in Chapter 2, crossover discrepancies have 
been the primary observable in the reduction of the orbit error in 
altimeter data from Geos3, Seasat and currently Geosat. Adjusted sea 
surface heights have been used to create global mean sea surfaces 
and estimates of the sea surface topography. These surfaces 
obviously contain all radial orbit errors that are not observable from 
crossovers. 

If we assume that all the signal in a crossover discrepancy is due 
to radial orbit errors then we can write 


Ar 0 = Ar(t!) - Ar(t 2 ) (5.106) 

where t! represents the time of the ascending crossing and t 2 the 
time of the descending crossing. For simplicity it is also assumed 
that the crossing occurs between arcs of the same integration period. 
Then the discrepancy using the first order periodic error will be 


Ar Q = 2Arp - 2Ar* - 2Arji 


(5.107) 
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where *v* refers to the variable part of the various error components. 
The part of the error that cancels out during the differencing is 


Ar° = Ar° - Ar 0 - Ar° + Aa x + Ar x 


(5.108) 


where *c’ refers to the mean geographic error. It is composed of a 
constant bias (Aa x - Ar 0 ), of a latitude dependent function (Arf - Arf) 
and the periodic function Arp which is both latitude and longitude 
dependent. 

The spatial signatures of both the mean and variable parts of the 
first order periodic error of gravitational origin are shown in Figures 
22 and 23 respectively. These maps are constructed using the same 
potential coefficient errors as before and plotted on a 5* grid with a 
50 cm contour interval. We can observe that both types of errors 
have quite an erratic variation with individual values ranging from 
about -3 meters to 2 meters. They both have the same regional 
characteristics although with different magnitudes. These 

characteristics follow the directions of the subsatellite tracks and are 
prominent in regions of rapid across track variations of either 
ascending or descending arcs (compare with Figures 20 and 21). The 
particular patterns of these errors may very well change in shape 
and or locations with a different choice of potential coefficient errors 
and of course with a different initial integration point. 

Obviously, in a crossover adjustment, Ar c cannot be observed and 
remains intact in the adjusted surface. Judging from Figure 23 this 
error can be substantial and can affect the computation of the 
stationary SST and geophysical investigations even on a local basis. 
This situation becomes worse when we consider the time dependent 
errors. The variable part of these errors is primarily a function of 
the crossing latitude and the crossing time differences (t,-t 2 ) while 
the part that cancels out is simply Ar T (t!) (T denotes all the time 
dependent errors and t, is assumed to be smaller than t 2 ). This type 
of error may very well create a high frequency across track 
signature in the mean radial error, because spatially neighbouring 
crossovers may belong to pairs of arcs that are widely separated in 
time, resulting into large variations in small spatial distances. 

The dependency of the mean error on each particular arc can 
create additional problems when crossovers between different arcs are 
formed. Then, each crossover discrepancy (considering for simplicity 
only Ar Gl ) is 


Ar o = 2Arp - Ar^ 


Ar? ij - 


(Ar? i + Ar? j ) 


(5.109) 
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Figure 22. Mean Radial Error of a Seasat Arc Using the GEM10B- 
GRIM3L1 Half Differences as Potential Coefficient Errors. 
C.I. - 50 cm. 
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Figure 23. Variable Radial Error of a Seasat Arc Using the GEM10B- 
GRIM3L1 Half Differences as Potential Coefficient Errors. 
C.I. = 50 cm. 



88 


with additional terms for the initial state errors and the time 
dependent errors. In (5.109) i and j denote the two different weekly 
arcs, and i j denotes the difference of the corresponding mean errors, 
which are generally not the same from arc to arc (see Table 5). The 
unobserved part of the error will be 


Ar c = Ar^ - min(Aro, Ar„) - minCAr? 1 , Ar?^) 


(5.110) 


with additional elements for initial state induced errors. Formation of 
several crossovers in a small geographic area between several arcs 
(there are twelve weekly arcs for Seasat) will create a surface of 
unobserved error that has high frequency signature in both the 
along and across track directions. The magnitude of this signature 
clearly depends on the differences of the mean errors in the various 
arcs. The incorporation of the time dependent errors makes the 
situation even worse. Further study needs to be made for particular 
patterns of this error where crossovers between different weekly arcs 
are considered. 

An additional error in the geometric adjustment of the crossover 
discrepancies arises from the fact that one arc of about 30-40 minutes 
of length is kept fixed so that to remove the singularities existing in 
the adjustment. Then the error in this arc (both the mean and 
variable parts) remains unaltered and propagates through all the 
crossing arcs, which in turn propagate the resulting error to all the 
parallel arcs in an unpredictable way. 

These types of problems are expected to have occured in all 
determinations of mean sea surfaces (e.g., Rapp, 1986). The problems 
are enhanced when data from both Geos3 and Seasat satellites are 
combined, since mean errors in the two types of orbits are expected 
to differ appreciably. The resulting geographic mean error has 

clearly a high frequency signature that can be easily detected in the 
detailed maps of the sea surface and the gravity field that have been 
produced (ibid, 1986). 

In terms of potential coefficient corrections it can be seen in the 
equations that all degrees and orders are partly unobservable from 
crossover discrepancies. More particularly the observability of the 
zonal coefficients is zero from the variable part of Ar p and is very 
small in the constant and 1 cy/rev parts since all the coefficients are 
lumped together. The observability of the potential coefficient 
corrections increases with increasing order. 


CHAPTER VI 


ANALYSIS OF GEOID UNDULATION ERRORS AND STATIONARY SST 


In Chapters 4 and 5 the complete analysis of the radial orbit 
error was made and formulas expressing the error in a Lagrangian 
form, in a Fourier series form and a geographic representation, were 
developed. Furthermore, the theory has been tested for its accuracy 
and numerical results were given. A similar analysis can be made for 
the geoid undulation errors and the stationary SST. Since both 
quantities are expressed as a spherical harmonic expansion on a 
sphere, this analysis is common for both of them. Some differences 
between the two quantities exist and will be pointed out. The 
spherical harmonic expansion of the geoid undulations is 

N = I ff*] I (C* cosmX + S, sinmAlP. (sin*) (6.1) 

i 7 £ ~ 2 ^ * * in— 0^*^ ^ rn ) * m 


where (4>,A,r) are the geocentric coordinates of the computation point, 
7 is the normal gravity on the reference ellipsoid^ P are the fully 
normalized associated Legendre functions and Cf m , are the 

normalized coefficients of the disturbing potential. The C$ m 
coefficients are related to the potential coefficients by 


(">♦ = r — p re 

Urn SLm 


( 6 . 2 ) 


with all Cf m ref being equal to zero except for the even degree zonals. 
The summation from 1-2 in (6.1) implies that the reference ellipsoid is 
a geocentric ellipsoid having the same mass as the earth. The geoid 
undulations can be computed using potential coefficients of some 
gravity field model that is determined either from satellite data or 
from terrestrial and altimetrically derived gravity anomalies. As it is 
stated in Chapter 3 a satellite derived gravity field provides only the 
long wavelength undulation signatures. On the contrary, use of a 
detailed (l*xl* or even finer grid) set of gravity anomalies can 
provide potential coefficients of high degree and order which in turn 
can be used to compute high frequency undulations. Specific high 
degree models that have been computed are the OSU81 up to degree 
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180 (Rapp, 1981) and OSU86F up to degree 360 (Rapp and Cruz, 1986). 
Use of any of those models can provide undulation estimates which 
are related to the true undulations by 


N = N 0 + AN C + AN° 


(6.3) 


where N 0 are the computed undulations, AN C are the commission errors 
and AN 0 are the omission errors. Estimates from the OSU86.F field 
which is complete to degree 360 (although with some degree of 
smoothing in degrees higher than 180) indicate that the undulations 
up to harmonic degree 36 have an RMS magnitude of 30.43 meters with 
an omission error of 1.58 meters and a commision error of 79 cm. 
Undulations up to harmonic degree 180, on the other hand, have an 
RMS magnitude of 30.47 meters, an omission error of 23 cm and a 
commission error of 1 meter. This indicates that about 95 percent of 
the undulation signal is at wavelengths corresponding to harmonic 
degrees smaller than 36, while more than 99 percent of the signal is 
contained in degrees smaller than 180. These commission and omission 
errors are approximate estimates and are based on the assumption 
that there is no signal above degree 360 and that the accuracy 
estimates of the field are reliable. In reality the omission errors are 
expected to be slightly higher because of the smoothing of the field 
above degree 180 and the omission errors of the field itself, which 
though on the average should be quite small. 

The omission errors are obviously expressed by equation (6.1) with 
the summation over £ being from i max + 1 to infinity. The commission 
errors can be written 


_ , . ^ ma x f B i - - 

AN = -*= X l (AC. cosmX + AS. sinmXlP. (sin*) 

ry A = 2 1 r > m =o 

(6.4) 

where ACi m and ASi m are the potential coefficient errors. Equation 
(6.4) can provide the commission errors at the groundtrack of a 
satellite if the geocentric coordinates of the subsatellite points are 
given. These coordinates are given in the Geophysical Data Records 
and are computed from the fully perturbed state vector of the 
satellite. 

In order to be able to express (6.4) in an inertial coordinate 
system with respect to Keplerian elements, a transformation of the 
surface spherical harmonics on the orbital plane has to be applied. 
Then 
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. it *max f a 1 ^ ^ i ^ 

AN* = -^ Z I I AC. I F. (i) 

rr £=a l r ) m=0 a =o ^ma p=o ^mp v 


i 


cos[(i-2p) (w+f) + m(Q-0) + C] 


(6.5) 


where C is given by (5.12) and i, «, f, Q, are the osculating elements 
of the satellite with groundtracks 4> and X. Note that (6.5) is given 
with respect to unnormalized coefficients for consistency with the 
developments made so far. In order to obtain an expression 
consistent with the ones for the radial orbit error (equation 4.62 or 
5.29), a circular orbit approximation and a linearization in the context 
of section 4.5.1 are necessary. Considering also that the magnitudes 
of the potential coefficient errors are small, a spherical earth 
approximation can be made. Then (6.5) reduces to 


AN° 



IF. (i)cos[ (X-2p) (M+w) + m(Q-0) 

p *mp 


+ c] 


( 6 . 6 ) 


where now i, w, M and fl take values on the reference orbit and can 
be computed from equations (4.36). A test on the effect of these 
approximations on the computation of AN C was made by using the half 
differences of GEM10B - GRIM3L1 as potential coefficient errors. 
First, equation (6.4) was used to compute &N C along the subsatellite 
points of the first six day Seasat arc at two minute intervals. This 
error is plotted in Figure 24. It has an RMS magnitude of 1.15 
meters with maximum values on the order of 9 meters. The 
commission error was also computed using (6.6). The RMS discrepancy 
with the "true" values is 10 cm with a maximum discrepancy of about 
70 cm. The discrepancies between the two sets occur because of the 
linearization that is used in «, M and Q, resulting in computing the 
undulation error not at the correct coordinates but at a location 
implied by the secularly changing «, M and (1. This level of 
agreement indicates that equation (6.6) cannot be used to model the 
commission undulation errors. These effects have to be modeled using 
equation (6.4). 

Although equation (6.6) cannot be used for modeling purposes, it 
is very useful in investigating the frequency content of the 
undulation errors. For this investigation equation (6.6) has to be 
rearranged through an analysis similar to the one in section 5.3. 
Setting k = i-2p and reordering the summations we obtain 


Imax imax ^max J 

4N = \-l Jo li . Jo 1C <o.. F «o^‘ <i)COS(l( l <n, lt + * k „„ + C) 


(6.7) 
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where f km and f kmo are defined in (5.8) and (5.13) respectively, ^ min 
is given by (5.61) and the summation in A is over degrees that have 
the same parity with k. After some trigonometric operations we 
obtain 


AN C 

= E fa cos ip At + b sinip At) 

^ km km km km > 



(6.8) 

where 

a 

km 

kpar fj 

= R | F 4l ^Ci) [1 
2 

[ ■e.^mlcosY' + 

l-AS. i "kmo 

Am 

(£;-] 

Am 

simb 

' kmo J 

(6.9) 

b 

km 

k par r i 

= R \ [1 

2 

Eic( m ) C ° S Amo ‘ 1 

*m 

AS im 

*m 

sinV' 

' r kmoJ 

(6.10) 


or, in terms of amplitudes and phases 


AN = £ A cos (V' At - ) 

km km km 


( 6 . 11 ) 


where 


A 


km 



( 6 . 12 ) 




km 



(6.13) 


Equations (6.8) - (6.10) provide the Fourier series representation of 
the undulation commission errors. The frequencies are again 
expressed in terms of cy/rev and the phases are computed based on 
the mean elements « 0 , M 0 , 0 o of the satellite arc. Equations (6.8) - 
(6.10) have a similarity with the corresponding equations (5.57) - 
(5.59) for the zero order periodic radial error. The actual spectral 
signatures of the two quantities though are expected to differ 
appreciably, due to the damping factor applied to the inclination 
functions in equation (5.58), causing the frequencies that are close to 
zero and 1 cy/rev to have large amplitudes, and higher frequencies 
to have small amplitudes, as shown in Figure 15. 
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The amplitude spectrum of the commission error is shown in 
Figure 25 up to 6 cy/rev. It is observed that this error has a more 
or less flat spectrum with the exception of the integer cy/rev 
frequencies that have no energy. It turns out that the higher 
frequencies do not die out as in Figure 15 but they have energies up 
to 36 cy/rev, which is the maximum frequency defined when using a 
field up to harmonic degree 36. The RMS errors by degree and order 
have also been computed and shown in Figures 26 and 27 
respectively. Both Figures show distinct differences from the 
corresponding Figures 16 and 18. Figure 26 shows that indeed the 
GEM10B and GRIM3L1 models have large coefficient discrepancies at 
degrees 6 to 12. 

As discussed in Chapter 3 the modeling of the undulation errors 
up to a high degree is not practical and does not yield an accurate 
solution for potential coefficient recovery since residual sea surface 
heights can only be obtained in the oceans. So the same gravity 

field that is used for the orbit integration is also used for the 
computation of the reference undulations. It has also been clear that 
such a field, being of low degree and order (l max = 36) yields a large 
omission error that coexists in the residual sea surface heights 
together with the radial orbit errors, the commission undulation 
errors and other effects. This omission error is on the order of 1.58 
meters (as computed using OSU86F). The frequency content of the 
omission error can be determined by equations of the type (6.8) - 
(6.10). These equations have been implemented using the OSU86F 
potential coefficients from degree 37 to 180. The maximum degree 180 
is chosen because of the practical difficulty in computing the 
inclination functions for degrees higher than 180 even with the 
extremely fast method derived by Goad (1987). It turns out that the 
omission errors have a substantial effect both in the lower part of 
the spectrum (k<36 cy/rev) as well as in the higher part (k>36 
cy/rev). The RMS magnitudes of both these effects are shown in 
Table 7. 

The higher frequency omission undulation errors of 1.20 meters 
can in principle be filtered out by applying a low pass filtering to 
the sea surface heights with a cutoff frequency of 36 cy/rev. This 
filtering is expected to have problems close to the ocean boundaries, 
where the profiles become discontinuous. The lower frequency effects 
of 1.04 meters are more serious because they cannot be separated 
from the commission errors since all degrees of a particular order are 
lumped together to generate one frequency. Furthermore additional 
frequencies lower than 36 cy/rev are generated from harmonics of 
higher orders. Both these lower frequency errors can be reduced by 
modeling the higher degree undulations errors using a high degree 
and order gravity field (i.e. 0SU86F) and removing these undulations 
from the residual sea surface heights. This removal will introduce 
commission errors due to the uncertainties of the high degree field. 
These errors have been computed based on the standard deviations of 
the OSU86F field up to degree 180 and are shown in Table 7. We can 
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Figure 26. RMS Geoid Undulation Error by Harmonic Degree and 
Cumulative, Using the GEM10B - GRIM3L1 Half Differences 
as Potential Coefficient Errors. 
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Figure 27. RMS Geoid Undulation Error by Harmonic Order Using the 
GEM10B - GRIM3L1 Half Differences as Potential Coefficient 
Errors. 
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see that the errors both at high and low parts of the spectrum are 
substantially reduced. A low pass filtering is still needed to remove 
the high frequency errors. The side effects due to the profile 
discontinuities are now expected to be much smaller. The error of 32 
cm at the low frequencies is expected to affect the estimated potential 
coefficient corrections and the stationary SST coefficients. A 
quantitative estimate of this effect though has not been obtained. 
Further study is needed in this matter. 


Table 7 

RMS Geoid Undulation Omission Errors for a Gravity Field up to 
l max = 36 as Computed from the OSU86F Field and Its Standard 

Deviations to i max = 180 


Frequency 

Signal 

Error 

k < 36 cy/rev 
k > 36 cy/rev 
Total 

1.04 m 

1.20 

1.58 

0 . 32 m 

0.48 

0.60 


A similar analysis with the preceding one for the undulation 
errors, can be made for the stationary SST. The stationary SST can 
be viewed as a function that is zero on land and takes specific values 
in the oceans. Furthermore it can be approximated as being a 
function on a sphere with the mean radius of the earth. Then the 
periodic part of the stationary SST can be expanded into surface 
spherical harmonics 

® £ _ _ _ 

f = R E I 1C, cosm\ + S 4 sinmMP. (sin«t>) (6.14) 

£- l m=o ' * m T * m T 


where Ci mT , S^ mT are the spherical harmonic coefficients of the 
stationary SST on a sphere of radius R. The zero degree term, if 
any, is not present in (6.14) because it is absorbed by the mean 
earth ellipsoid that is used as a reference for the geoid undulations. 
This is consistent with the definition of the oceanic geoid which 
requires that the average SST be zero as sampled globally in oceanic 
regions. A more detailed discussion on the definition of the oceanic 
geoid and how the mean earth ellipsoid can be defined in the 
presence of the stationary SST is contained in Engelis (1985). 
Contrary to (6.1) equation (6.14) contains a first degree term which 
arises from the depression of the sea surface, relative to the geoid, 
in the southern hemisphere. This depression is caused by the strong 
circumpolar currents. 
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The transformation of (6.14) into a Fourier series form follows the 
same steps and has the same problems with the one leading to 
equations (6.8) - (6.10). But now i=l is also included giving rise to a 
1 cy/rev signature, in contrast to undulation errors that do not have 
any 1 cy/rev effects. Modeling of the stationary SST in the context 
of the combined altimeter problem has the same problems regarding 
the omission errors. Solution up to degree 36 introduces aliasing 
effects into frequencies up to 36 cy/rev arising from higher degree 
harmonics. The higher frequency signature due to those harmonics 
will be filtered out together with the higher frequency undulation 
errors. 



CHAPTER VII 


COVARIANCE ESTIMATION OF RADIAL DISTANCES AND GEOID 
UNDULATIONS BASED ON A GEOPOTENTIAL COVARIANCE 


7.1. Introduction 

The numerical example that was worked out in Chapters 5 and 6 
has provided error estimates for both the radial distances (for a 
Seasat orbit) and undulations. These error estimates, although they 
are realistic and may indeed be a possible choice of errors that might 
occur, they are only based on an arbitrary set of numbers 
representing potential coefficient errors. In reality, errors in 
existing gravity field models are unknown. The only available 
information in some of these models is the estimates of their variances 
and, for even fewer of them, their error correlations. So the only 
information that one can infer about the errors of the radial distances 
and undulations is of statistical nature. Using the covariance matrix 
of a gravity field model one can compute variances and error 
correlations of the above quantities and of their frequency 
components by 


l = GI G T 

C S 


(7.1) 


where I cs is the covariance matrix of potential coefficients and G is 
the Jacobian that can be computed in a straightforward manner, since 
all the models developed so far are linear with respect to potential 
coefficients. 

Since the orbital specifications that have been used so far in this 
study are the ones of the Seasat orbits as computed from the PSSS4 
gravity field, it was decided to also use its covariance matrix to 
retain a consistency throughout the study. This will also provide 
some accuracy estimates for computed orbits that have already been 
used for the computation of mean sea surfaces (e.g. Marsh et al., 
1986a). As discussed in Chapter 2, the PGSS4 gravity field is given 
complete to degree and order 36 and has been determined from the 
combination of the PGSS3 gravity field with a set of Seasat altimeter 
data. Such a tailoring of the field to the Seasat characteristics led to 
the computation of orbits with a reported RMS error of 70 cm (Lerch 
et al., 1982) which was based on the computed RMS crossover 
discrepancy of two 6 day Seasat arcs. 
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For the computation of reliable accuracy estimates and error 
correlations of both the orbit and the geoid an appropriate scaling of 
the covariance matrix has to be made. Such scaling is common 
practice for the most recent GEM models (Lerch et al., 1985) but has 
not been applied on the PGSS4 matrix. For the purposes of this 
study the PGSS4 covariance is scaled up by a factor of 10. This 
scaling being uniform in all degrees and orders is suboptimum and 
may result into erroneous accuracy estimates. But since the primary 
purpose of this study is to demonstrate the ability of analytic 
techniques to model and eventually determine the radial orbit error, 
the very accurate scaling of the covariance matrix is not that critical 
at this stage. 

In the next sections the accuracies of the radial distances and 
undulations as well as their error correlations are going to be 
computed both with respect to time and in their geographical 
representation. For the temporal representation the Fourier series 
expression is going to be used since it is more efficient and 
economical than the Lagrangian form. The Fourier series expression 
will also enable us to compute accuracies of individual frequencies 
and their correlations as well as the accuracies of the radial distances 
by degree and order. It should be noted that secular and second 
order effects are not going to be computed. The geographic 
representation on the other hand will provide the accuracy estimates 
of ascending and descending arcs as well as statistical estimates of 
the parts of the orbit error that are observed or non observed in 
crossover discrepancies. 


7.2. Covariance Propagation Using the Fourier Series Approach 

For the efficient computation of (7.1) the Fourier series 
expression has first to be converted in a matrix form. Equation (5.19) 
can be written as follows 


Ar 


km 


cosii' At sinf At 

km km 


km 


km I 


where 


c 


km 


S 


km 



cos^ sintf 

r km r km 

-sin-tf cos ip 

km km 


D 

E 


km 

km 


(7.2) 


(7.3) 
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From (5.24) and (5.25) D k(n » E km can be written 
D = l fP .AC. - Q, .AS. ) 

E = I fQ .AC. + P, .AS. } 

km A m in ' kmA *m kmA 

where 


(7.4) 

(7.5) 


P kmA 


Q kmA 


A . 
kmi 

-B . 

km^ 

B J 

km^ 

A * 

km^ 


&—m even 


&— m odd 
&— m even 


^*~m odd 


(7.6) 


(7.7) 


& . = raax(lk!-2, m, 2) (7.8) 


and A km i, B km i are given by (5.16) and (5.17). In a matrix form we 
have 


D 

km 


P 


km<£ 


mi n 


P - -Q • • . . -Q 

km^ m * km-G • km^ 

max. min r 


AC, 


mi n n 


AC 


^max m 


AS, 


mi n n 


AS, 


(7.9) 


with a similar expression for E km . 


In a more compact notation 
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D 

P 

- Q 

AC 

m 

km 

km 

km 

AS_ 

m 




*c m 

— 

Q - P 

*S m 


km km 


(7.10) 


(7.11) 


where P km , Q km , AC m , AS m are row and column vectors of length 4 max - 
£ min + 1* The generalization for all frequencies km can be obtained 
in a straight-forward manner. Then the matrix equation for all the 
coefficients C km , S km can be written in the following compact form 


c 


km 

- 

S 


km 



cosV 1 sinafr 

r km ~ km 

-sin i> cosy 

~ km r km 





m 


m 


(7.12) 


where C km and S km are vectors of length (2* max + 5)(i max + 1). Since 
the energy in the orbit error is negligible in frequencies higher than 
6 cy/rev each of those vectors has a length of 481. The coefficients 
C km and S k(n are stored orderwise from 0 to 36 and for k starting 
from -6 to 6. The quantities cosy^ m and sin^ km represent diagonal 
matrices of dimension 481x481. Each of the diagonal elements contain 
the cosine or sine of the phase that corresponds to that particular 
row. and AS£ m are vectors containing the potential coefficient 

corrections. The vector AC has [ (^+l)(4+2)/2-3] (or 700) elements 
since no zero and first degree terms are present, while the vector 
ASi m has a length of [i(i+l)/2-l] (or 665) elements since in addition 
no zonal terms are present. All coefficients are arranged orderwise 
so that to be consistent with the structure of the PGSS4 covariance 
matrix. Finally the matrices Pg m and Qg ro referring to AC£ m have a 
dimension of 481x700 and the matrices Pjj m and Qg m referring to AS^ m 
have a dimension of 481x665. All four of these matrices have mostly 
zero elements except for (* max - * mjn + 1) elements per row that are 
non-zero and have values as given by (7.6), (7.7) and (7.9). Equation 

(7.12) readily provides the Jacobian (or the design matrix) of the 
Fourier coefficients of the radial orbit error. 

To compute the covariance matrix of the Fourier coefficients we 
first compute the covariance matrix of the coefficients D km and E km 
for one particular combination of indices km and nj. Equation (7.1) 
then gives 
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(7.13) 


where the geopotential covariance of orders m and j is split into the 
covariance matrices Ec» Es> Ecs> and Esc* Then 


, = P Z P + Q E Q -PE Q - Q I P 

km,nj km C nj km S nj km cs nj km SC nj 


= CT E.Q . + P 5 Z P s + Q c z P s + P s z Q c 

’ km , n j c n J km S nj km CS nj km SC nj 


r D E 


= P IQ. - Q I P - Q I Q + P Z P 

km,nj krn ^ nj km S nj km SC nj km CS nj 


(7.14) 

(7.15) 

(7.16) 


Since most of the elements in these matrix multiplications are zero it 
is more efficient to express these equations in a summation form over 
the non-zero elements. Then we obtain 


r n - I I fP *P . . ffr + Q a Q ..<Tc 

km,nj £ . I ^m, i j n J^ ^m, i j 

- P <fpc - Q A P <r r c 1 

km* nji CS i m ,ij km* nji CS ij J i m J 


(7.17) 


r F = Z I TQ «Q ffr + P .P <Tc 

km,nj £ . L km * *m,ij kmt "j 1 *m, i j 

+ Q A P <Tpc + P *Q (fpc 1 

km* nji CS * m ,ij km* nji CS ij,£ m ] 


(7.18) 


r 0 e - Z I [P. »Q . ,<r 

km,nj £ . I 


L, I * .v< up Q .P <f e 

. [ km* nji km* nji 5, ^ (llj \ j 


+ P km* P nji <Tcs i m> j j °km* Q n j i * c s j j > £ J (7.19) 


The above equations can compute the variances of the coefficients D k 
and E km when km=nj, covariances between coefficients of different 


jsr 3 
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or different n, or as computed above, covariances between coefficients 
of any frequency. From equations (7.17) - (7.19) we can compute the 
covariances of the Fourier coefficients as follows 


= cost/) cosy ff D + siny sin^ .ov 

'km,nj n J km,nj * tm n J km,nj 


+ ( cosy, siny . + siny cos y .)<r n p 

km nj km nj km,nj 


(7.20) 


s. = sin^ sin^ ni <r 0 _ + cos* cos* <r E 


km, n j 


n J km,nj km n J km,nj 


- (cosY' sinY . + sinV' cosY' ,)<r n c 

km nj km nj km,nj 


(7.21) 


r r c = -cosy, sinf ,ff n + siny cosy ,<r c 

km,nj * cm n j km,nj n j km,nj 


+ (cosy cosy - sin* sin* )<r n|: 

~ km nj r km n j DE km , n j 


(7.22) 


Computation of all the covariances (7.20) - (7.22) generates the 
covariance matrix of the Fourier coefficients. This matrix is a 

symmetric one and can be split into three matrices namely I c , £ s , and 
Ecs> where X c and I s are also symmetric. Then the variance of a 
coefficient of frequency km and its covariances with all the other 
elements can be found at the k row where 


(<W+k) (k-x+k+l) 

k = <..x(2k.. x +l)(k m . x *k) j + » + 1 (7-23) 


Its covariance with a particular nj Fourier coefficient is found in the 
column that has the following index with respect to the diagonal of 
the R row 


index = n - k + 1 


(7.24) 


where n is computed from (7.23) with n, j as arguments. The above 
is valid for either E c or E s . Since £ cs is a square matrix the 
elements are simply located by the indices h and it. This indexing is 
very important in the data management aspect of the computation and 
usage of this matrix. 
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The determination of the variances of the Fourier coefficients at 
zero and 1 cy/rev frequencies as well as their correlations with all 
the other frequencies is now a trivial matter, since their Jacobian can 
be readily obtained from equations (5.31), (5.33), (5.34) and (4.54) - 
(4.61). The basic difference between this Jacobian and the one that 
led to equations (7.20) — (7.22) is that all its row elements are non 
zero since all the potential coefficient terms give rise to these 
frequencies. Furthermore this Jacobian is clearly arc dependent in 
contrast to the first one that is arc independent and representative 
of the accuracies of the orbit in general. 

A similar reasoning exists for the resonant and second order 
induced errors. Accuracies and error correlations due to zonal 
coefficients can be easily computed using the same formulation and 
considering the covariances for the zonal coefficients, while all the 
second order effects are readily available since they are just linear 
functions of the constant bias Aa 0 and the resonant terms. 

Having the covariance matrix of the Fourier coefficients the 
variance of the radial distances at time At can be computed by using 
the following generalization of (7.2) 


Ar 


I lcosV' k(n At sinV , k(n Atl 



(7.25) 


where now the cosine and sine are diagonal matrices of dimension 
481x481 with their elements computed for the frequencies km and a 
particular At. Obviously C km and S km are row vectors containing all 
the Fourier coefficients and I is a unit row vector. Implementation of 
(7.1) can then provide the variance of the radial distance at At. With 
the same reasoning we can obtain the covariance matrix of all or a 
subset of radial distances within an arc or among different arcs. 

The computation of the statistics of the uncertainties of the 
radial distance is based on equations (5.37), (5.38). As a matter of 
fact it turns out that the equations describing the RMS estimates of 
the uncertainties are identical to (5.39) - (5.41), applied for the 

variances. Then the RMS radial orbit uncertainty is 



(7.26) 


The RMS uncertainty by order becomes 
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Ar„ 



(7.27) 


while the RMS uncertainty by degree becomes 


t 






(7.28) 


By breaking (7.28) even further we can obtain the RMS uncertainties 
for any degree and order. 

The variances of the Fourier coefficients up to 6 cy/rev have 
been computed by implementing (7.20) - (7.22) and using the scaled 
PGSS4 covariance matrix up to degree 36. The corresponding 

spectrum of uncertainties is shown in Figure 28, while the RMS radial 
orbit uncertainty is computed to be 63 cm. From Figure 28 wo 

observe that the uncertainty at 1 cy/rev is very small. All the 
uncertainties are practically at frequencies below 2 cy/rev. It is also 
seen that the frequencies very close to 1 cy/rev have very small 
uncertainties which are due to the tailoring of the field with the 
Seasat altimeter data. Recognizing that this tailoring affected 
primarily the harmonics above degree 20 the following test was made 
to understand what this effect has been. The variances of the 
Fourier coefficients have been recomputed using the scaled PGSS4 
covariance up to degeee 20 (Solution II) and the covariance between 
degrees 21 and 36 (Solution III). A final variance computation used 
the covariances between the coefficients of degrees lower than 20 
with those of degrees greater than 20 (Solution IV). It is obvious 
that the sum of these three computations is equal to the first 
computation that used the full covariance up to degree 36 (Solution 
I). Comparison of the individual variance components between the 
four solutions for all frequencies, indicates that for most of the 
frequencies solutions III and IV give no significant contributions. In 
other words the use of the scaled PGSS4 covariance up to degree 20 
gives almost identical variances as with the use of the covariance 
complete to degree 36. An exception to the above is for four 
frequencies that are very close to 1 cy/rev. These frequencies 
contain the first harmonic order and the resonant orders 14 and 15. 
They are shown in Table 8. The variance components of solutions II, 
III, and IV for all the four frequencies are such that the total 
variances of solution I are very close to zero. These components are 
shown in Table 8. 



flDIRL 


Table 8 
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As it can be seen from Table 8 the incorporation of the Seasat 
altimeter data has introduced significant correlations between 
coefficients of orders 1, 14 and 15 so that to produce lumped 

coefficients at frequencies very close to 1 cy/rev, that give an almost 
perfect fit to the orbit. 

The RMS uncertainties of the arc dependent constant bias and 1 
cy/rev effects have also been computed and are shown in Table 9 for 
all the six day arcs of Seasat. The total RMS uncertainty for each 
arc has also been computed, taking into account the arc independent 
RMS uncertainty of 63 cm, and is shown in Table 9. 


Table 9 

Constant, 1 cy/rev and Total RMS Radial Distance Standard Deviations 
for all Seasat Arcs Based on the Scaled PGSS4 Covariance. 


Seasat Arc 

Constant 

1 cy/rev 

Total 

1 

0.39 m 

0.70 m 

1.02 m 

2 

0.28 

0.64 

0.94 

3 

0.38 

0.67 

0.99 

4 

0.38 

0.70 

1.01 

5 

0.39 

0.67 

0.99 

6 

0.48 

0.72 

1.07 

7 

0.58 

0.77 

1.15 

8 

1.07 

1.08 

1.64 

9 

0.38 

0.67 

0.99 

10 

0.39 

0.71 

1.03 

11 

0.66 

0.85 

1.25 

12 

0.37 

0.69 

1.01 


In all the arcs the 1 cy/rev error is the dominant error. It has 
an RMS uncertainty generally on the order of 70 cm with the 
exception of the eighth and eleventh arcs. The constant bias is 
around 40 cm with higher values for the above two arcs. No 

explanation for these higher uncertainties has been found. The total 

RMS uncertainty is on the order of 1 meter, again with the exception 
of the eighth and eleventh arcs. The accuracy estimates in Table 9 
agree quite well with the 0.70 m accuracy estimate that was given by 

Lerch et al. (1982), if we consider that the latter estimate is based on 

the computed RMS crossover discrepancy for two 6 day arcs and that 
the part of the orbit error that is not observed by crossovers has 
almost the same magnitude with the part that is observed. So if we 
assume that the non observed part also provides an RMS uncertainty 
of 0.70 m, then the total RMS accuracy is 0.99 m which is consistent 
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with the estimates of Table 9. This agreement indicates that the 
scaling of the PGSS4 covariance by a factor of 10 could be realistic 
although, it is definitely a simplified one. As mentioned at the 
beginning of this Chapter, these accuracy estimates do not include 
any long wavelength zonal effects or second order effects. This 
should not cause any problems because these effects are expected to 
have been properly removed from the Seasat orbits. Additional errors 
arising from initial state vector errors, air drag, solar radiation 
pressure and other error sources are not expected to significantly 
change these estimates. 

Additional statistics that have been computed include the arc 
independent RMS accuracies by degree, by order, and by degree and 
order. These are shown in Figures 29 and 30 and in Table 10, 
respectively. Figure 30 shows again the small uncertainties in orders 
1, 14 and 15 which result from the use of the altimeter data. From 

Figure 29 we observe that the odd degrees give the largest 

contribution to the orbit error, together with degrees 14, 15, and 16. 
Note that the overall uncertainty from Figure 30 is much smaller than 
in Figure 29 because of the significant correlations between 

coefficients of different degrees and same order which are not 
accounted for in Figure 29. Basically the same type of information 
can be obtained from Table 10 where we can observe the large 
uncertainties of individual odd degree coefficients and primarily the 
ones of harmonic orders 1, 14, and 15. Again the correlations 
between these coefficients are not accounted for. So Figure 29 and 
Table 10 cannot really be used by themselves for an accuracy 

assessment since they provide a quite distorted estimate of the orbital 
accuracies. 

In order to see what is the behavior of radial accuracies with 
respect to time, a computation was made of the variances of the radial 
distances for the first Seasat arc, at time intervals of 5 minutes, and 
using the scaled PGSS4 covariance up to degree 36. Correlations 

between radial distances for the first three revolutions have also 
been computed. To ease the computational burden, only frequencies 
up to 2 cy/rev have been used. The full covariance matrix of the 
Fourier coefficients up to that frequency as well as the variances of 
the constant bias and 1 cy/rev effects and their covariances with all 
the other frequencies have been computed for this purpose. The 
accuracies for the first 6 day Seasat arc are plotted in Figure 31 
where we can observe the strong 1 cy/rev effect and daily and 
semidaily terms. Note that at the beginning of the arc the 

uncertainty is very close to zero and it accumulates rapidly at the 1 

m level. The small error at 3 days occurs because the satellite 
passes through the same geographic region after that period of time 

and so the error is similar to the error at the initial minutes of the 

arc. This is another indication of the geographic correlation of the 
orbit error. 
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Figure 29. RMS of Seasat Radial Distance Standard Deviations by 
Harmonic Degree Based on the Scaled PGSS4 Covariance. 



HARMONIC ORDER 

Figure 30. RMS of Seasat Radial Distance Standard Deviations by 
Harmonic Order Based on the Scaled PGSS4 Covariance. 
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Table 10 


of Seasat Radial Distance Standard Deviations for each Coefficient 
Degree and Order Based on the Scaled PGSS4 Covariance. 

Units are in Centimeters. 
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Figure 31. Seasat Radial Distance Standard Deviations Based on the 
Scaled PGSS4 Covariance. 



Figure 32. Error Correlations Between Seasat Radial Distances Based 
on the Scaled PGSS4 Covariance. 
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Figure 32 shows a representative correlation curve for time lags 
up to 300 minutes (approximately 3 revolutions). The correlation 
obtains a maximum negative value of about 0.40 after half a revolution 
and a local maximum positive value of about 0.70 after one revolution. 
This pattern continues with smaller amplitudes in the correlations, to 
become zero between radial distances being 1.5 days apart or 
equivalently differing by 180* in longitude. After that, the 
correlations start increasing again to become very close to 1 for a lag 
of 3 days. The actual magnitude of the correlations depends on the 
particular time that we consider as initial point for the correlation 
computation, but is not expected to vary appreciably. 

The covariance propagation of geoid undulations can be made in 
an identical way to the one for the radial distances. All the 
equations developed so far in this section can be used for this 
purpose with the following changes in (7.6), (7.7), and (7.8) 


P 


kmfi 


i-n 


i— m odd 


(7.29) 


Q 


km.£ 


RF. fi-k 

*> m 


n 


fi — m odd 


(7.30) 


and fi min as given by (5.61). The computations have to be made for 
all frequencies up to 36 cy/rev since the undulation spectrum has 
significant energies at the higher frequencies also. In a similar 
analysis the covariances of the errors in the residual sea surface 
heights can be obtained. In this study only the variances of the 
undulations and residual sea surface heights have been computed up 
to 6 cy/rev. The corresponding spectra are shown in Figures 33 and 
34. Note again the distinctly different behavior of the radial distance 
(Figure 28) and undulation spectra. The residual sea surface height 
spectrum (without the sea surface topography signal) has amplified 
uncertainties at frequencies below the 1 cy/rev and smaller ones at 
frequencies between 1 and 2 cy/rev. This happens because the error 
correlations between the radial distances and undulations are positive 
in the first frequency band and negative in the second frequency 
band. 


7.3. Covariance Propagation Using the Geographic Representation 

Variances and covariances of radial distances and geoid 
undulations can also be obtained using equations (5.86), (5.87), (5.91), 
and (6.4). The covariance propagation in terms of geographic 
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coordinates has the advantage of providing the spatial behavior of 
the accuracies in a satellite orbit as well as the geoid. Consistent 
with the computations in section 5.4 the accuracy estimation is to be 
done on a geographic grid. Then interpolation to the actual point of 
altimeter observation can provide the accuracy of the orbit at that 
point. As described in section 5.3.1, the first order radial orbit error 
of gravitational origin can be written 


Ar = Ar c * Ar v (7.31) 

with 

Ar° = Arj; - Arg (7.32) 

Ar v = Ar^ - Arg (7.33) 

where 

Arg = Ar 0 + Arj (7.34) 

Arg = Arg (7.35) 


and c, v denote the mean and variable parts of the error 
respectively. Then the variances of the above quantities become 


__ 2 _ c 2 . „ y2 ± 

<r r - <r r + <r r 


(7.36) 

_c 2 - _C 2 4. _.c 2 

- ff p + *0 

- 2ff po 

(7.37) 

2 - 2 . -V 2 

2<r P° 

(7.38) 

ff cv = ff cv + ,cv 

- 

(7.39) 


To compute all the variances and covariances in equations (7.36) - 
(7.39) requires some tedious algebraic manipulations. Here only the 
final expressions for each of these quantities are going to be given. 
To simplify the notation we denote 



+ B, 


&mp 


c *. P ) Q ir -<*> 


tmp 


(7.40) 
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A L = *( A * mp + B W + (7 - 41) 

(7.42) 

(7.43) 

(7.44) 

(7.45) 

(7.46) 

(7.47) 


A* = 

( B i„ 

+ c ): 

' . )cosmX 0 + (b. + C. 

sim; 0 l £n) £ m 

sin4>] 

siniJ 

sinmX 0 

(7.48) 

B. = 

K 

+ °im 

sin *]sil«X < , + [b? + c“ 

SiniJ ° K &m 

sin<t>1 

siniJ 

cosmX 0 

(7.49) 

O 

to 

3 

II 

(i- 

sin 2 ^ 

sin 2 i 

1 ^ [ 2 v 1 v 

C. cosmX 0 + C* smmX 0 

1 l. xi n 



(7.50) 

O 

to 

3 

II 

[>- 

sin 2 4 l l 

sin 2 iJ 

S f 2v 1 v ] 

C. sinraX 0 + C. cosmX 0 

v -Km -Km J 



(7.51) 


B, = X A. Q. (4> 0 ) 

J>m p *mp *mp 


B 


m 


1 1 v, r «. p ( *» ) 


,1c 

'Am 

,2c 

,2v 

'im 


= I Wt (♦„) 

p *">P 

= i X W* (♦„) 

p * mp 

= X w* (♦„) 

p *">P 


i X W 4 __ (♦„) 

p 


^mp 


then 


. 2 ^max ^max r ^max ^max r r 

= X X X X ,,f2-tf. . 1A? A .cosmXcosjX<r c . 

m =o 1=m i=i" 1 im ij J c *m,ij 


^max ^max c c 

+ X I (2 - 6 t 1A. A. . sinmXsinjX<r s 
j=m i = i" l jJ £m i j J s .t m , 


i J 


^max A max c c 

+ 2 I I A. A. .cosmXsinjX<r cs 

j = l i=i^ * m ’J ^m, i j 


(7.52) 
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where 


t = max (2,m) 
i = max (2,j) 


i if j=m 
j if j>m 


1 if i=& and j=m 
0 otherwise 


Similarly 


ffp = II n f 2-6. )A* A v . sinmXsin jX<r c 

P 1 1 «m, i ji tm l j J«m, i j 


+ I I (2 -6. lA, A v cosmXcosjX<r s 

' ij *m,ij 


- 21 Z A* A v sinmXcosjX<r cs 

*® ' j £m, i j 


*7 = l l[l E j C 0 S ®^ sin j^"*‘( 1— i j) A^A^ jSinmXeos jX j < 

- £ A v .sinmXcosjX+(l-6^ A^A^cosmXsinjXJi 

- E X(A^ A v cosmXcosjX-A^ A c ,sinmXsinjXj<r cs ^ . j 

<ro = l l[l 1(2-6. 1 (A. A <r c . +B . B <r s . ) 

*■ l i j J ^ ij c ^m, i j ij s ^m, i j J 


+ 2X X A B. ,<r cs . ( .] 

*m i j *m, i j J 


(7.53) 


(7.54) 

c ^cn , i j 

s ^m, i j 

(7.55) 


(7.56) 
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r* = I l\l 1(2-6. ) fC. C <r c , +D, D <r s . ) 

° L l ^mij s ^m,ijJ 

+ 2££ C i. D J 
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(7.57) 


(7.58) 


(7.59) 


+ 1 1 ( A L Il iJ siI “ x - c « m A (j cosjx )'’'=*j„,,j] 


(7.60) 
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1 J ( A I„ D , jSinjX)...^ , j] 


(7.61) 


-t; = n[n (*I.A, J *i»ft + (l-«,. jlJ )A 1 ,XjS l nj x )**l., t i 
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(7.62) 


In equations (7.54) - (7.62) all the summations are over the indices 
indicated in (7.52). For the uniformity of equations the quantities 
referring to the zonal S coefficients 


5 *o, 


1 j 


= 0 



= 0 


have also been included. The above equations provide the variance 
of the radial distance at a point (4>, X) as well as the variances 
representing the mean and variable radial errors. The same type of 
equations can be used to compute spatial covariances for radial 
distances at different geographic locations and/or between different 
integration arcs, by assigning the appropriate values to the latitude 
dependent functions, the cosines and sines of longitudes and the arc 
dependent functions. 

The undulation variance can be computed using an equation 
similar to (7.52) but replacing the functions A c by the Legendre 
functions. Then the residual sea surface height variance can be 
computed by the combination of the above equations. 

A correct implementation of the equations with the complete scaled 
PGSS4 covariance, requires the computations to be made on at least a 
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1* grid. This was assessed by computing the variances on a small 
area and then interpolating and comparing to the results of the 

previous section. A global implementation though on a 1* grid 
becomes extremely time consuming and has not been made. It has 
been estimated that about 20-25 CPU hours on an IBM 3081D 

mainframe are required for this task even for well optimized software. 
For demonstration purposes computations have been made on a 5* 

grid using the scaled PGSS4 covariance up to degree 20. It turns 
out that even for that maximum degree the grid of 5* is not adequate 
because it does not provide accurate interpolated values, particularly 
in the northern and southern latitudes. Therefore the accuracy maps 
plotted in Figures 35-39 are only approximations to the correct 
accuracies. Note in all Figures that the plotted uncertainties are 
much higher than the uncertainties computed in the previous section 
where the full geopotential covariance up to degree 36 was used. 

Also in Figure 39 that shows the undulation accuracies we can see 
that the oceanic regions are represented much more accurately than 
the continents. Both features are yet another indication of the 
tailoring of the PGSS4 field with the Seasat altimeter data. 



Figure 35. Standard Deviations of Radial Distances of Seasat 
Ascending Arcs Based on the Scaled PGSS4 Covariance to 
Harmonic Degree 20. C.I. = 10 cm. 
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Figure 36. Standard Deviations of Radial Distances of Seasat 
Descending Arcs Based on the Scaled PGSS4 Covariance to 
Harmonic Degee 20. C.I. = 10 cm. 
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Figure 37. Standard Deviations Reflecting the Geographic Mean Error 
of Seasat Arcs Based on the Scaled PGSS4 Covariance to 
Harmonic Degree 20. C.I. = 10 cm. 
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Figure 38. Standard Deviations Reflecting the Geographic Variable 
Error of Radial Distances of Seast Arcs Based on the 
Scaled PGSS4 Covariance to a Harmonic Degree 20. 
C.I. = 10 cm. 


Figure 39. Geoid Undulation Standard Deviations Based on the Scaled 
PGSS4 Covariance to Degree 20. C.I. = 10 cm. 











CHAPTER VIII 


IMPROVEMENT OF THE ORBIT AND DETERMINATION OF THE 

STATIONARY SST 


8.1. Introduction 

For the solution of the combined altimeter problem the parameters 
Asj, ip and ip T of equations (3.13) or (3.15) have to be optimally 
determined. As discussed in the previous Chapters the computation 
of the partials of (3.13) can be made by numerical integration for the 
radial distance related partials, and using equations (6.4) and (6.14) 
for the geoid undulations and stationary SST. Choosing to work with 
the analytical approach, it has been seen that the Lagrangian theory 
is quite accurate in modeling the radial orbit errors, and gives a 
considerably better insight to its properties. It has also been seen 
that the geographic representation of the radial error is not very 
accurate when it is to be implemented on a regular grid and it 

becomes extremely expensive for computation on a point by point 
basis. Furthermore, it has turned out that the Fourier series form is 
the most practical and economical to use without any loss of accuracy. 
On the contrary, undulation errors and stationary SST cannot be 
accurately modeled by a Fourier series expansion and equations (6.4) 
and (6.14) have to be implemented. 

The detailed observation equations can be easily obtained by 

using the expressions for the Jacobian of the errors as derived in 
Chapter 7, complemented by the partials corresponding to the 

resonant and second order terms as well as the initial state vector 

errors. These equations are not going to be given here since this 
would be a repetition. Instead, problems associated with the 
estimability of the coefficients are going to be discussed and ways to 
improve it will be proposed and applied. 


8.2. Problems Involved in the Estimation 

From the Fourier series expressions of radial orbit error, geoid 
undulation error and stationary SST it is obvious that the potential 
coefficient corrections and the SST coefficients have several types of 
dependencies. More specifically the potential coefficients of all the 
degrees and of a particular order are linearly combined to create 
lumped coefficients that give rise to specific frequencies, each 
separated by 1 cy/rev and modulated by m cy/day. For a gravity 
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field up to degree and order 36 there are 77 such lumped coefficients 
for each order, providing a total of 2849 lumped coefficients which 
under ideal conditions (no other errors or SST signal, optimum 
sampling, etc.) can in principle provide a solution for the 1365 

potential coefficients contained in the field. Even in such a case 
there are problems that prohibit an efficient estimation. The most 
serious problem arises from the considerable correlations between the 
potential coefficients of same order and the correlations that exist 
between the frequencies themselves. Furthermore due to the parity 
contraints in the undulation expansion (and partly in the orbit error 
expansion) the degrees of freedom of the system decrease 

substantially. Finally, as it has been seen in Chapters 5 and 6, 

frequencies corresponding to zonal terms do not really have any 
significant amplitude, resulting into low estirnability of these 

coefficients. Coefficients of higher degree and order have a greater 
degree of freedom, since they are involved in linear combinations of 
fewer terms, and so they are expected to be better estimated. The 
same types of dependencies are valid for the stationary SST 
coefficients. 

The estirnability problem is further complicated by the coexistence 
of the three types of unknowns that give rise to the same 
frequencies. More specifically the initial state vector errors are fully 
correlated with the 1 cy/rev errors and the second order errors of 
gravitational origin. On the other hand, the stationary SST has the 
same spectrum as the combined spectrum of orbit and undulation 
errors at least in most of the frequencies. This can be seen from 
Figures 33 and 34. Some separability can be obtained only for the 
low frequencies (lower than 1 cy/rev). Furthermore the first degree 
zonal terms of the stationary SST are fully correlated with the 1 
cy/rev gravitational and initial state vector errors. 

Higher degree (< max > 36) geoid undulation effects create 
additional problems. As it has been seen in Chapter 6 even if we 
assume that a perfect filtering of the higher frequency undulations 
can be obtained (with the implementation of higher degree fields and 
low pass filtering) there are still higher degree undulation effects 
that are aliased into low frequencies. 

Finally the altimeter observations are not sampled globally but 
only in oceanic regions. This gives rise to additional errors in the 
coefficients to be estimated. These errors are due to leakage effects 
resulting from the discontinuities of the observables over land. 


8.3. Conditioning of the System and Solution 

To overcome the problems discussed in the previous section, or at 
least to improve the estirnability of the unknowns, a certain 
conditioning has to be applied to the normal equations resulting from 
(3.13) or (3.15). Such conditioning can be obtained by introducing 
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prior information for the parameters that can be of two types. 
Strong prior information which translates to constraining some of the 
parameters that are considered to be very well known, or weak prior 
information, which consists of providing information about the 
statistical properties of the parameters. 

For the potential coefficient corrections strong information can be 
obtained by constraining the low degree zonals which are considered 
to be known quite accurately. In such a case the long period effects 
of zonal origin, particularly in the eccentricity, have to be properly 
removed. Statistical information on the other hand can be provided 
by the covariance matrix of the gravity field that is used for the 
orbit generation. This matrix, being practically the major means to 
separate the unknowns, has to be a reliable and well scaled matrix, 
accurately depicting the apriori correlations of the coefficients. A 
similar covariance matrix can also be obtained for the initial state 
vector of each arc from the orbital adjustment. The stationary SST 
coefficients can be conditioned by use of oceanographic information. 
A rather mild conditioning can be obtained by using the degree 
variances resulting from the spherical harmonic expansion of existing 
SST models (e.g. Levitus, 1982). Such a spherical harmonic expansion 
has been carried out by Engelis (1983, 1985). 

A further conditioning that can be applied to the solution is to 
use the crossover discrepancies as observables, so that to provide an 
additional type of observation equations. The mathematical model for 
a crossover discrepancy is 


p(t 2 ) - p(t x ) = r(t 2 ) - h(t 2 ) - r(ti) + h(t x ) = 0 (8.1) 


where tj and t 2 are the times of crossing between two arcs. The 
convention followed in this differencing is that t 2 is greater than t,. 
Such a convention allows for the easier computer implementation of 
the observation equations. Changing the time attributes with the 
corresponding indices and using (3.9) - (3.14) we obtain 


dr d§ 3r d§ t 

ds 2 ds T dSj dSj 

( 8 . 2 ) 


Asi 


dr 2 ds. 
3 5 3n 


dr t ds ] 
35 3n 


Ap = 0 


where it has been assumed that there are no time variations of the 
sea surface or any other time varying effects. This assumption is 
consistent with the development of equations (3.13) and (3.15). 
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Equation (8.2) models the differences of the radial orbit error at 
two epochs and is independent of any geoid undulation and stationary 
SST signature. So crossover discrepancies can in principle be used 
by themselves to solve for potential coefficient corrections and initial 
state errors, in order to correct the computed orbit. Due to the 
differencing though, the observability of all the lumped coefficients is 
reduced. This observability becomes null when the interval t 2 -ti of 
the differenced sea surface heights is an integer multiple of the 
period of the frequency generated by each lumped coefficient. For 
example, the 1 cy/rev error is completely unobservable at the 
northernmost and southernmost latitudes where crossover 
discrepancies have time differences very close to multiples of one 
revolution. This error component is perfectly observable at the 
equator. This results in a reduced sensitivity of equation (8.2) to 
most of the parameters. In particular, zonal coefficients cannot be 
estimated from crossover discrepancies. 

Similar to sea surface heights, crossover discrepancies also have 
dependencies between their parameters. Again initial state vector 
errors are fully correlated with the 1 cy/rev and second order 
gravitational errors. Furthermore, potential coefficient corrections 
within the same lumped coefficients are highly correlated. So, 
conditioning with the covariance of the field is needed to improve the 
estimability of the parameters. But even with the conditioning the 

estimated parameters are such that, although they do minimize the 
crossover discrepancies, they do not necessarily minimize the radial 
orbit error also. As can be seen in the next section, there is indeed 
a substantial reduction of the radial orbit error but it is not as 
effective as when sea surface heights are also used. Furthermore no 
simultaneous determination of the stationary SST coefficients can be 
obtained. So crossover discrepancies can be efficiently used only 
when combined with sea surface heights. In such a combination, 
crossover discrepancies can provide an excellent conditioning of the 
system since they do not contain any undulations or SST signature. 

The observation equations for the residual sea surface heights 
and the crossover discrepancies have the following form respectively 


V t - A p Ap T Aj 


Ap 

Ap T 

As t 


Ah 


V 2 = | Ap Aj 


Ap 

As, 


+ (h 2 - h t ) 


(8.3) 


(8.4) 


or in a more compact notation 
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V = AX + L (8.5) 

A very important aspect in forming the observation equations is 
the sampling of the observations. According to the sampling theorem, 
in order to be able to resolve frequencies up to 36 cy/rev, the 

residual sea surface heights have to be sampled every one minute for 
a time period of approximately six days (for the case of Seasat). The 
sampling of the crossover discrepancies is more, complicated. The 

reason is that sampling has to be made over time differences and not 

the time itself. Definitely the consideration of all the crossover 

discrepancies is suboptimum since crossover formations are more 
dense in the northern and southern latitudes than in the equatorial 
regions and so they "observe" the gravity field in a non-homogeneous 
way. This irregular sampling becomes even worse when crossovers in 
repeat arcs are considered, so appropriate downweighting has to be 
applied. No quantitative investigation has been made on the sampling 
of the crossovers and more study is felt to be necessary. 

The solution of the equations (8.5) with prior information is a 
typical weighted parameters least squares solution and is given by 


X = - (A t PA + P x )“ 1 A t PL 


( 8 . 6 ) 


where P is the weight matrix of the observations and P x is the 
inverse of the prior covariance for the unknowns. In a typical 
solution for a six day Seasat arc there are about 2500 unknowns to 
be estimated, using over 10000 observations. The combination of the 
two data types can create singularities if both the sea surface 
heights h! and h 2 at a crossing location and the corresponding 
crossover discrepancy ha-hi are used as observables since these 
three quantities are linearly dependent. When either h t or h 2 and 
h 2 -hi are used, their error correlation has to be considered. 

An additional outcome of the solution is the a posteriori 
covariance matrix of the parameters which can yield error estimates 
for the orbit, the sea surface, the stationary SST and the long 
wavelength geoid, as well as their error correlations. 


8.4. A Simulated Solution 

As discussed in the previous section, the effectiveness of the 
solution heavily depends on the consistency of the prior information 
with the unknown parameters, and on the types of observables to be 
used and their sampling, while additional factors like a priori 
accuracy of the orbit, precision of the altimeter observations and 
others are also quite important. So simulated solutions taking into 
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others are also quite important. So simulated solutions taking into 
account all these factors will be made to assess the optimum 
conditions for the setup of the solution and its sensitivity to these 
factors. 

From the discussions made so far, it is obvious that the solution 
of the combined altimeter problem is computationally extremely 
intensive. The computer resources that are required both in terms of 
computing time and memory are typical of a supercomputer 
environment. For the purposes of software development and 
validation of the method it was felt that such intensive computations 
would be an unnecessary burden. So a computationally more relaxed 
setup was considered in which, simulated residual sea surface heights 
and crossover discrepancies up to harmonic degree 10 for a 3 day 
Seasat arc were used. 

For the generation of the simulated "unknown" potential 
coefficient corrections the scaled PGSS4 covariance matrix up to 
degree 10 was used to create multivariate normally distributed random 
deviates. Using these deviates, radial orbit errors of gravitational 

origin (including second order effects) and undulation errors were 

computed every 3 minutes for the first 3 days of the first Seasat arc. 
No initial state vector errors were considered. The harmonic 

coefficients of the Levitus SST up to harmonic degree 10 were used 
to compute the stationary SST at the same intervals. Then, residual 
sea surface heights were generated by adding a 10 cm white noise to 
account for the altimeter noise. For the generation of crossover 

discrepancies, the radial orbit errors (with 10 cm white noise) were 
computed and differenced at the arc crossing times. The algorithm 
that is used to determine the crossing times as well as the geocentric 
latitudes and longitudes of the crossings is described in Appendix E. 
A similar algorithm is given in Shum (1983). The total number of 
observations for this 3 day Seasat arc were 984 residual sea surface 
heights and 1178 crossover discrepancies. For the weight matrix of 
the parameters, the inverses of the scaled PGSS4 covariance to degree 
10 and of the degree variances of the Levitus SST field were used. 

Using the above observations, three solutions were made. The 

first solution used both the residual sea surface heights and 
crossover discrepancies, the second solution used only the residual 
sea surface heights and the third only the crossover discrepancies. 
The recovered coefficients, for all three solutions, were then compared 
to the initial "true" coefficients. The RMS percentage discrepancies 
by degree and order for each of the three solutions (identified as 
solutions 1, 2, 3) are shown in Tables 11 and 12, together with the 
"true" RMS magnitudes in terms of geoid undulation corrections and 
stationary SST. 
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Table 11 

RMS Magnitudes and Percentage Discrepancies by Degree Using the 

Scaled PGSS4 Covariance 


Degree 

Undulation Correction 

Stationary SST 

Magn. 

Sol 1 

Sol 2 

Sol 3 

Magn. 

Sol 1 

Sol 2 

1 





18 cm 

59% 

79% 

2 

4 cm 

43% 

41% 

27% 

22 

6 

12 

3 

7 

26 

24 

51 

9 

28 

36 

4 

6 

31 

45 

58 

7 

30 

33 

5 

9 

52 

53 

77 

7 

67 

68 

6 

12 

19 

30 

30 

14 

29 

21 

7 

12 

46 

53 

64 

10 

42 

46 

8 

13 

28 

31 

36 

7 

50 

49 

9 

18 

24 

24 

30 

5 

90 

83 

10 

13 

i 

22 

37 

43 

3 

119 

155 


Table 12 

RMS Magnitudes and Percentage Discrepancies by Order Using the 

Scaled PGSS4 Covariance 


Order 

Undu! 

.ation Correction 


Stationary 

SST 

Magn. 

Sol 1 

Sol 2 

Sol 3 

Magn. 

Sol 1 

Sol 2 

0 

2 cm 

97% 

100% 

124% 

30 

cm 

23% 

23% 

1 

13 

48 

48 

55 

16 


61 

77 

2 

17 

29 

38 

55 

8 


46 

67 

3 

12 

45 

40 

63 

7 


77 

71 

4 

10 

17 

32 

46 

3 


98 

98 

5 

8 

41 

32 

38 

4 


110 

87 

6 

13 

8 

38 

19 

2 


68 

138 

7 

13 

2 

7 

5 

2 


77 

92 

8 

7 

4 

12 

12 

1 


76 

107 

9 

4 

5 

28 

4 

1 


110 

67 

10 

4 

3 

6 

6 

1 


110 

118 
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In an overall evaluation, the estimation of the coefficients has 
been satisfactory. The recovery of the potential coefficient 
corrections is generally better than the coefficients of the stationary 
SST. This is due primarily to the heavier weighting of the former 
ones which allows for a better decorrelation and less influence by the 
altimeter noise. The estimability of the potential coefficient 
corrections by degree is almost the same for all the degrees 
(discrepancies are between 20-50%). Orderwise though there is a 
great variation, starting from practically no estimability for the zonal 
terms to perfect recovery for the higher order terms. The SST 
coefficients being only mildly constrained are more affected by the 
noise. This can be seen in particular in the higher degrees and 

orders where the energy of SST is small and the RMS discrepancies 
become higher. At low degrees, and particularly the second degree, 
the SST coefficients are recovered very well with the exception of the 
first degree and order terms, which are fully correlated with the 1 
cy/rev radial error. 

The recovered potential coefficient corrections and SST 
coefficients have been used to compute radial orbit errors, commission 
undulation errors, crossover discrepancies, stationary SST and 
residual sea surface heights. All these quantities were then compared 
to the "true" ones over the oceanic regions. The RMS "true" 
magnitudes and their corresponding discrepancies are shown in Table 
13 for all three solutions. 


Table 13 

RMS True Magnitudes and Discrepancies Using the Scaled PGSS4 

Covariance 



True 

Magnitude 

Discrepancy 



Sol 3 

MB 

111 cm 

11 cm 

15 cm 

19 cm 

1 

36 

10 

12 

14 


50 

15 

18 

— 


114 

4 

5 

— 

Ar D 

106 

3 

5 

3 

^as c 

88 

11 

14 

19 

^ r desc 

97 

11 

14 

19 


There are many conclusions that can be drawn from Table 13. 
Again it is obvious that the first solution that uses both types of 
observables provides much better estimates. In terms of the radial 
orbit error the recovery is at the 90% level. The recovery of 
undulation corrections and stationary SST is on the order of 70%. 
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The residual sea surface heights are perfectly recovered, although 
individual discrepancies of Ar, AN, and { are much higher. This is 
due to the correlations existing in these quantities, particularly 
between Ar and AN. The crossover discrepancies are effectively 

minimized although the individual errors in the ascending and 
descending arcs at the crossing locations have the same RMS 
magnitude as the RMS radial error. 

The solution using only residual sea surface heights performs 
slightly worse than the first solution and is better than the third one 
that uses only crossover discrepancies as observables. In this last 
solution it can be seen that although the crossover discrepancies are 
perfectly recovered, the individual radial errors are worse than in 
solution 1 by a factor of 2. The . undulation correction recovery has a 
similar degradation. The most important disadvantage of using 
crossover discrepancies is that the stationary SST cannot be 
determined simultaneously. A suboptimum solution can be obtained 
using methods described in Engelis (1985). Then the adjusted values 
of the sea surface and geoid can provide a stationary SST that would 
have an RMS discrepancy from the "true" values of about 25 cm. 

The original error magnitudes shown in Table 13 are rather 
pessimistic considering that they only reflect errors up to degree 10. 
Nevertheless the estimation reduces these errors quite drastically. In 
order to see what is the degree of improvement when a more accurate 
gravity field is used the simulated solution was repeated using the 
unsealed PGSS4 covariance. The RMS magnitudes and discrepancies 
for all quantities are shown in Table 14. 


Table 14 

RMS True Magnitudes and Discrepancies Using the Unsealed PGSS4 

Covariance 



True 

Magnitude 


|| 

Sol 1 

Sol 2 

.. 

Sol 3 

■S^Hj 

35 cm 

6 cm 

7 cm 

6 cm 

1 

10 

5 

5 

5 

B 

50 

8 

9 

- 


55 

4 

5 

- 

Ar □ 

33 

3 

4 

3 

Ar aso 

28 

6 

6 

6 

^ r desc 

31 

6 

6 

6 













131 


From Table 14 it can be seen that in this case all solutions give 
practically the same results although the first one is still superior. 
Note that use of crossover discrepancies now gives better results 
than the residual sea surface heights. The percentage of recovery of 
the errors, as compared to the previous solution, is smaller and is 
subject to a limit imposed by the altimeter noise. On the contrary, 
the stationary SST, in the presence of smaller errors, is determined 
much better .than in the previous solution. 

These simulated solutions can show the apparent requirement for 
a very accurate gravity field in order to efficiently compute the 
stationary SST from satellite altimetry. They also show that even 
very accurate altimeter data like the Topex/Poseidon altimeter data 
can further be improved. Finally they show that altimetric data from 
any mission even of moderate accuracies can be substantially 
improved and that quite reliable SST estimates can be derived. A 
necessary requirement is the existence of a well scaled and reliable 
geopotential covariance. To test the importance of this requirement 
several other solutions were made. In one of these solutions the 
PGSS4 variances were only used. In a second solution the PGSS4 
covariance was slightly altered while in other ones potential 
coefficient errors that were completely inconsistent with the PGSS4 
covariance were used. Other solutions did not use any prior 

information for the stationary SST. In all solutions the results were 
worse than the ones contained in Tables 13 and 14 and in many of 
them the results were unacceptable. 

The "true" quantities (or a priori errors) and the discrepancies 
after the adjustment (a posteriori errors) as well as the Levitus and 
the estimated SST have also been computed on a geographic basis on 
a regular 5” grid. Additionally the a priori and a posteriori 
accuracies based on the scaled PGSS4 covariance and the estimated 
covariance, have been computed on the same grid. All these estimates 
are shown in Figures 40-61. The first conclusion is that the radial 
errors both for the ascending and descending arcs have been 
substantially reduced. Figures 40 and 41 show the substantial 
reduction of the error of the ascending arcs in the Pacific Ocean 
where errors of the order of 3.5 meters have been reduced to 16 cm. 
On the contrary the errors in the South Atlantic Ocean have only 
been reduced by a factor of two. A similar improvement exists for 
the descending arcs (Figures 44 and 45) where errors of about 2 
meters in the South Atlantic Ocean have been reduced to about 20 cm. 
Figures 48 and 49 show the substantial reduction of the mean 
geographic error, while Figures 52 and 53 show that the variable 
error has been completely eliminated. Comparing Figures 41, 45, 49 
and 53 it becomes obvious that the residual radial orbit errors are 
common to both ascending and descending arcs. The apriori and 
aposteriori standard deviation maps in Figures 42, 43, 46, 47, 50, 51, 
54 and 55 are absolutely consistent in magnitude and shape with the 
maps of the corresponding errors, indicating that they can be used to 
obtain reliable estimates about the errors themselves. 
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The apriori and aposteriori geoid undulation errors are shown in 
Figures 56 and 57. Most of the error signature existing in Figure 56 
has been eliminated. The only significant residual errors are in the 
Pacific Ocean with an amplitude of 20 cm and alternating sign. The 
standard deviation maps, shown in Figures 58 and 59, also indicate 
the improvement achieved from the adjustment. The determination of 
the stationary SST as shown in Figures 60 and 61 has also been 
successful. It can be seen that all the features of the apriori SST 
have been recovered both in shape and magnitude. 
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Figure 40. A Priori Radial Error of Seasat Ascending Arcs Based on 
Simulated Potential Coefficient Errors to Harmonic Degree 
10. C.I. = 20 cm. 



Figure 41. A Posteriori Radial Error of Seasat Ascending Arcs Based 
on the Discrepancies Between Simulated and Recovered 
Potential Coefficient Errors to Harmonic Degree 10. C.I. = 
2 cm. 
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Figure 42. A Priori Standard Deviations of Radial Distances of Seasat 
Ascending Arcs Based on the Scaled PGSS4 Covariance to 
Harmonic Degree 10 . C.L = 20 cm. 



Figure 43. A Posteriori Standard Deviations of Radial Distances of 
Seasat Ascending Arcs Based on the Estimated 
Geopotential Covariance to Harmonic Degree 10. C.I. = 2 
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Figure 44. A Priori Radial Error of Seasat Descending Arcs Based on 
Simulated Potential Coefficient Errors to Harmonic Degree 
10. C.I. = 20 cm. 
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Figure 45. A Posteriori Radial Error of Seasat Descending Arcs Based 
on the Discrepancies Betwen Simulated and Recovered 
Potential Coefficient Errors to Harmonic Degree 10. 
C.I. = 2 cm. 
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Figure 46. A Priori Standard Deviations of Radial Distances of Seasat 
Descending Arcs Based on the Scaled PGSS4 Covariance to 
Harmonic Degree 10. C.I. = 20 cm. 
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Figure 47. A Posteriori Standard Deviations of Radial Distances of 
Seasat Descending Arcs Based on the Recovered 
Geopotential Covariance to Harmonic Degree 10. C.I. = 2 









137 


0 30 60 90 120 150 180 210 210 270 300 330 0 



-70 >r> rr i= r ■ f *ri- 7o 


0 30 60 90 120 150 180 210 240 270 300 330 0 


Figure 48. A Priori Mean Geographic Radial Error of a Seasat Arc 
Based on Simulated Potential Coefficient Errors to 
Harmonic Degree 10. C.I. = 20 cm. 
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Figure 49. A Posteriori Mean Geographic Radial Error of a Seasat Arc 
Based on the Discrepancies Between Simulated and 
Recovered Potential Coefficient Errors to Harmonic Degree 
10. C.I. = 2 cm. 
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Figure 50. A Prion Standard Deviations Reflecting the Mean Geo- 
graphic Radial Error of a Seasat Arc Based on the Scaled 
PGSS4 Covariance to harmonic Degree 10. C.I. = 20 cm. 
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Figure 52. A Priori Variable Geographic Radial Error of a Seasat Arc 
Based on Simulated Potential Coefficient Errors to 
Harmonic Degree 10. C.I. = 20 cm. 
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Figure 53. A Posteriori Variable Geographic Radial Error of a Soasat 
Arc Based on the Discrepancies Between Simulated and 
Recovered Potential Coefficient Errors to Harmonic Degree 
10. C.I. = 2 cm. 
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Figure 54. A Priori Standard Deviations Reflecting the Variable 
Geographic Radial Error of a Seasat Arc Based on the 
Scaled PGSS4 Covariance to Harmonic Degree 10. C.I. = 20 
cm. 



Figure 55. A Posteriori Standard Deviations Reflecting the Variable 
Geographic Error of a Seasat Arc Based on the Estimated 
Geopotential Covariance to Harmonic Degree 10. C.I. = 
2 cm. 
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Figure 58. A Priori Standard Deviations of Geoid Undulations Based 
on the Scaled PGSS4 Covariance to Harmonic Degree 10. 
C.I. = 5 cm. 



Figure 59. A Posteriori Standard Deviations of Geoid Undulations 
Based on the Estimated . Geopotential Covariance to 
Harmonic Degree 10. C.I. = 2 cm. 








CHAPTER IX 


CONCLUSIONS, RECOMMENDATIONS 


In this study a solution to the combined altimeter problem (as 
defined in Chapter 2) has been proposed. This solution differs from 
earlier solutions, that were of geometric nature, since it considers the 
dynamic properties of the satellite orbit. Sea surface heights and 
crossover discrepancies are used in a minimum variance solution to 
simultaneously determine radial orbit corrections, the sea surface and 
the stationary SST. 

An analytical approach for modeling the radial orbit error has 
been adopted. Using the Lagrangian theory a linearized expression 
for the error, supplemented by second order effects, has been 
derived. A limited evaluation of this expression has shown that the 
analytic approach can be accurate enough to model the radial orbit 
error. Alternative forms have been derived, that provide a better 
insight into the properties of the error. It turns out that the 
Fourier series approach is much more convenient to use and equally 
accurate to the Lagrangian form. The geographic representation is 
valuable in showing the spatial characteristics of the radial error but 
is of inferior accuracy due to approximations that are necessary for 
its implementation. 

It has been seen that the radial orbit error is of long 
wavelength nature. The most dominant part is at a frequency of 1 
cy/rev. A constant bias and errors in frequencies very close to 1 
cy/rev have also a large magnitude. Additional error of 1 and 2 
cy/rev with a time dependent amplitude can also arise due to 
resonant conditions in the satellite orbit and due to second order 
effects. These time dependent errors can grow up to significant 
magnitudes if not properly removed during the orbital adjustment. In 
the space domain the radial error shows a considerable variation, that 
is systematic and depends on orbit specifications and the epoch of 
the orbit integration. It also turns out that one part of the error is 
common to both ascending and descending arcs while the other part 
is of equal magnitude but with a different sign for the two arcs. 
This fact is the primary reason why empirical crossover adjustments 
have failed to really remove the radial orbit error from the altimeter 
data. 

In order to use the sea surface heights as observables the geoid 
undulations and stationary SST must be modeled. It has been seen 


144 


145 


that modeling of these quantities with respect to time in an inertial 
coordinate system is not very accurate due to approximations that are 
involved. This formulation, although not accurate enough, can show 
that the spectra of undulation error and stationary SST are identical 
and that they both differ significantly from the one of the radial 
orbit error. These spectra have significant amplitudes in all the 
frequencies. Therefore a solution to the combined altimeter problem 
requires that the higher frequency signal be removed. This can be 
obtained by using a high degree gravity field and a low pass 
filtering. 

For the solution of the problem, prior information is required to 
reduce the dependencies that exist among the parameters and improve 
the stability of the system. The covariance matrix of the gravity 
field, used for the orbit generation, and oceanographic information 
can be used for this purpose. A limited simulated solution has shown 
that the radial orbit error and the stationary SST can be effectively 
recovered. Statistical estimates from the solution are consistent with 
the simulated values. 

As discussed in the introduction, the motivation and the initial 
stage of this research work have been based largely on investigations 
and conclusions obtained by Colombo (1984). A similar analysis of the 
radial orbit error has been made by Wagner (1985) and Rosborough 
(1986). Wagner has derived formulas that are identical to the ones 
contained in Section 5.3 and he has obtained conclusions that are 
very similar to the ones of the present analysis. Rosborough has 
made a complete analysis of the radial orbit error both in the Fourier 
series form and its geographic representation. He has given 
numerical estimates of the error for the Topex/Poseidon satellite using 
the differences of GEM10B and GRIM3B gravity fields as potential 
coefficient errors. Furthermore he has derived a formulation that 

leads to the computation of statistical estimates of the error. He has 
used the GEML2 geopotential covariance complete up to degree 20 to 
provide numerical results for the Topex/Poseidon radial uncertainties. 

Both Wagner’s and Rosborough’s formulations are equivalent to 
the one presented here, although differences in the derivations are 
apparent. Neither of the two investigations has properly identified 
the 1 cy/rev error which is the largest error component in the 
determination of the radial distance to a satellite. This error 
component has been correctly determined in the present analysis, as 
seen in Section 4.7. An additional contribution of this study in the 
analysis of altimeter data, is the systematic development of the 
observation equations that lead to an optimum determination of the 
stationary SST and the identification of all the potential problems that 
are involved in such a determination. A similar estimation technique 
has been proposed by Wagner (1986) who has used approximate 
observation equations to model the residual sea surfce heights and 
only mild prior information to condition the unknowns. His simulated 
results indicate a lower estimability than the one presented in the 
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present simulation study. 

Several problems have not been addressed in this simulation 
study. The first and most important one involves a more complete 
modeling of the radial orbit error, using a complete gravity field up 
to degree 36. Such a solution was not made because of large 
computing requirements. Based on results from the solution up to 
harmonic degree 10 though, it is expected that the higher degree and 
order potential coefficient corrections can be very well determined. 
An exception may be the resonant order coefficients and so more 
testing needs to be made. 

Another problem that has to be addressed is how much the lower 
frequency high degree undulations contaminate the solution. This 
assessment can be obtained by simulating sea surface heights with 
the geoid undulations and stationary SST being modeled up to higher 
degree and order (e.g. 180). Then, using a high degree gravity field 
and low pass filtering the high frequency signal can be removed. 
The residual sea surface heights can then be used for the solution. 
This solution can be compared to a solution that does not contain any 
high degree signal. Conclusions on the contamination of the solution 
can then be obtained. 

Another simulation test that has to be made is related to the 
optimum sampling of the crossover discrepancies. If all the crossover 
discrepancies are used, then there is an oversampling at the northern 
and southern regions relative to the equatorial regions. This 
irregular sampling results into a non uniform observability of the 
gravity field. An optimum sampling can be obtained by applying the 
sampling theory for time differences. Such an application has not 
been made and more investigation is necessary. 

Even with all the above problems addressed, a solution of the 
altimeter problem as proposed in this study is oversimplified, 
primarily because of the incomplete modeling of the radial orbit error. 
As discussed in Chapter 2, radial orbit errors of smaller magnitude 
also arise from air drag and solar radiation pressure modeling errors. 
These errors affect all the Keplerian elements and primarily the mean 
anomaly. According to Kaula (1966) these effects are linearly 
dependent with time and so they closely resemble resonant and 
second order errors of gravitational and initial state vector origin. 
Further complications in the solution are expected from errors in the 
modeling of tides and the other sea surface variations. Finally, the 
modeling of the stationary SST and the prior information used are 
also simplified. It is expected that if more comprehensive models of 
the stationary SST, based on dynamic equations of motion of the 
oceans, are used, then a better estimability of the ocean parameters 
can be obtained. 

Even the most comprehensive solution taking into account all 
these problems cannot recover high frequency stationary SST because 



147 


of the filtering that is required. These high frequency effects can 
be obtained in local high energetic areas by combining the recovered 
sea surface estimates with estimates of local detailed geoids. These 
detailed geoids can be computed using gravimetric methods, if 
accurate and dense gravity observations exist. 

Based on the analysis made so far, it is expected that the 
analytical methods will be able to efficiently process real altimeter 
data. It is anticipated that the GEMT1 field, along with improved 
models for tides and surface forces and more accurate station 
coordinates will be used to compute new orbits for Seasat. These 
orbits will be merged with altimeter observations to provide more 
accurate sea surface heights. From initial tests, using the GEMT1 
field, the crossover discrepancies are on the order of 70 cm (Marsh et 
al. 1986b) implying an accuracy of about 70 cm for the sea surface 
heights (assuming that the mean geographic error has the same order 
of magnitude as the variable geographic error). This data set, being 
more accurate than the PGSS4 set will have the additional advantage 
of being supplemented by a well scaled and accurate geopotential 
covariance, as the GEMT1 covariance is expected to be. Using this 
dataset a solution can be made, for a six day arc to obtain adjusted 
sea surface heights, an improved long wavelength geoid and long 
wavelength SST. The author anticipates a posteriori accuracies on 
the order of 10-15 cm or even better for all the above quantities. 
The solution has to be repeated for all six day arcs of Seasat. 
Comparison of the SST estimates between solutions can provide 
variability of the sea surface, while averaging of the individual 
solutions can provide stationary SST over the lifetime of the satellite 
and an even more reliable long wavelength geoid. 

Application of the method to the Geosat data is expected to be 
rather problematic since the Geosat orbits during the Extended Repeat 
Mission of the satellite have accuracies on the order of 4 meters 

(Cheney et al., 1987). These orbits are computed using the GEM10 

gravity field. With such poor orbits and with no geopotential 
covariance available, the only possible analysis that can be made is a 
crossover analysis along the lines of Chapter 8. The resulting sea 
surface heights are expected to have inferior accuracies as compared 
to the ones for Seasat. Furthermore, no stationary SST can be 
obtained, while the determination of SST variations might be affected 
by the across track high frequency variations of the orbit errors. A 
considerably improved analysis of the Geosat data can be obtained if 
the GEMT1 field is used for the orbit determination of Geosat. 

For geophysical applications, that require a very detailed 

knowledge of the sea surface, combination of observations from 
different satellites can be made. The procedure described in this 
study can be used to independently adjust observations from 

different satellites and derive the ocean effects at different epochs. 
These effects, being generally different, have to be removed from the 
adjusted surfaces to provide estimates of the geoid. Since it is not 
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expected that the a posteriori accuracies of each individual solution 
will be the same, the geoid of inferior accuracy can be fitted to the 
one of superior accuracy to provide a more detailed and consistent 
surface. Such a fit can be made by a regular geometric crossover 
adjustment. 

The activities in the area of satellite altimetry have been 
substantial during the last ten years and are expected to culminate 
with the launch of the Topex/Poseidon satellite. In view of the 
unprecedented accuracies that are expected from that mission, both in 
terms of orbit determination and altimeter observations, it is expected 
that comprehensive algorithms will be required for the optimum 
reduction and analysis of the data. In this context, the author hopes 
that this study has contributed in a positive way to the development 
of such algorithms. 
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APPENDIX A 


SATELLITE STATE VECTOR TRANSFORMATIONS 

Transformation between inertial rectangular coordinates, 
Keplerian elements (elliptic coordinates) and earth fixed spherical 
coordinates are considered. 

A. Conversion from elliptic to rectangular coordinates. 

Al. Use Kepler’s equation 


M = E - esinE 


(A. 1) 


to compute the eccentric anomaly E 
A2. Compute q 



qi’ 


a(cosE-e) 
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A3. Compute position vector ? 
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A4. Compute q 
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A5. Compute velocity vector ? 

r*i 


r = 




R 3 (-Q)M-i)Rs(-«)q 


liJ 

B. Conversion from rectangular to elliptic coordinates. 

Bl. Compute radial distance r and tangential velocity v 

r = Irl = (x 2 + y 2 + z 2 )^ 

v = Irl = (x 2 + y 2 + £ 2 )^ 

B2. Compute angular momentum vector fi 


h = ? x r 
hi = yz - yz 
h 2 = zx - zx 
h 3 = xy - xy 

h = Ihl = (hf + h 2 + h?) y * 

B3. Compute radial velocity r 

- (»• - 3 s 

B4. Compute semimajor axis a and eccentricity e 
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B5. Compute eccentric anomaly E, true anomaly f, and mean 
anomaly M. 


0 a-r 

cosE = 

ae 

(A. 12) 
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(A. 13) 

* * ^‘(S3) 

(A. 14) 

* - ta "-(^-e iDE l 

(A. 15) 

M ~ E - esinE 
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B6. Compute right ascension of ascending node Q and inclination i 


0 = tarrl bfe) 


(A. 17) 
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B7. Compute argument of perigee u 
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(A. 19) 


« + f = tan 


ipii 


(A. 20) 


u - (w+f) " f 


(A. 21) 



158 


C. Conversion from elliptic to spherical coordinates. 

Cl. Compute the eccentric anomaly E and true anomaly f according 
to (A.l) and (A. 15) 

C2. Compute radial distance r 

r = a(l-ecosE) (A. 22) 

C3. Compute geocentric latitude 4' 

4« = sin -1 [sinisin(cj+f)] (A. 23) 


4> is in the same quadrant as w+f if -tt/2 * <j+f * tt/2, otherwise in 
the same quadrant as rr - (o+f). 

C.4 Compute longitude X 


X = ■in-'( a>,i y t) ) +0-9 (A. 24) 


where X - 0 + 0 is in the same quadrant as sgn(7r/2 - i)(w+f) for 
X * tt/2 • When i = tt/2 then X - Q + 0 = 0 for cj+f < 7 t and 
X — Q + 0 - 7T for w+f - rr. (A.23) and (A. 24) are the equations of 
the groundtrack. 



APPENDIX B 


PERTURBATIONS DUE TO J 2 

In this Appendix perturbations in the elements a, e, M and their 
time derivatives due to J 2 are going to be computed. From equations 
(4.5) we have 
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The forcing function F is replaced by V 2G 
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(B . 3) 


V 


20 




G 2p q(e) cos [ (2-2p+q)M + (2-2p)«] 

(B.4) 


Then the derivatives of V 20 with respect to M, co, e, a can be 
computed as follows 
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Expanding the summation in p and q and considering that 


F200G20-I = ^ 202^221 (B. 6 ) 
F200G200 = F202G220 (B. 7 ) 
F201G21-1 = F201G211 (B. 8 ) 
F200CI2OI = F 2 02 G 22~1 (B. 9 ) 


we obtain 
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9 V 

the final expression for ** is 

oM 
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Expanding the summation in p, q we obtain 
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g *1 2 [ g] [2F 2 ooG 20 o cos(2M+2cj) 

+ 2F200G20— i cos(M+ 2 cj) + 2 F 2 oo^ 2 oi cos( 3 M~ 2 cj) 
+ 2F 2 o 1 ^ 2 1 — 1 cosM + F 201 G 210 ] 

( 1-e 2 )“ 3 / 2 and GJ 10 = 3e(l-e 2 )~ 5 / 2 


_ it 
a 


J 2 ft) [(| sin J i - |)e(l-e 2 ) 


2 \ ~ 5 / 2 


15 3 

+ — ^ esin 2 i cos(2M+2cj) + g sin 2 i cos(M+2a>) 

- sin 2 i cos(3M+2«) + sin 2 i - |]cosm] 


(B. 


29'. 


(B. 23) 


(B . 24) 


3^ 


I aJ 


I F 2 op G 2 Pq cos [ (2-2p+q)M + (2-2p)«] (B.25) 

p> q 
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Then 


aVj 

3a 


= |r J > ft)’ [(S sin=i - i) (l-e’r’/i 


3 3 

4 sin 2 icos(2M+2cj) + ^ esin 2 icos(M+2<a) 


- esin 2 icos(3M+2«) + sin 2 i-^j ecosM] 


(B.26) 


Using (B.17), (B.20) , (B.24) and (B.26) in (B.l), (B.2), (B.3) 
and setting 1 -e 2 ~ 1 for small eccentricities, we obtain 

a 20 = - 2na J 2 [^] sin 2 isin(2M+2«) - ^ esin 2 isin(M+2w) 

+ ^ esin 2 isin(3M+2«) - sin 2 i-^j esinNlj (B.27) 


e 20 = - nJ 2 ^ sin 2 isin(M+2cj) + ^ sin 2 isin(3M+2cj) 

- sin 2 i-jjj sinNlj (B.28) 

M 20 = nJ 2 [~(| sin 2 i-|j + ^ sin 2 icos(2M+2u) 

+ ^ sin 2 icos(M+2w) - sin 2 icos(3M+2cj) 

+ ^ (| sin 2 i-|] cosm] (B.29) 

where = n has been used in all three equations. In (B.29) all 

the terms that were multiplied by the eccentricity have been 
neglected. 
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Integration of (B.27), (B.28) and (B.29) for the periodic part of 
M 2 o> on the reference orbit gives 

a 20 = 2aJ 2 sin 2 icos(2M+2cj) - ^ esin 2 icos(M+2<j) 


+ S3 
24 

esin 2 icos(3M+2w) - sin 2 i -^jcosNlj 

(B.30) 

e 2 o = J 2 

(^a) [I s i n2 i cos (M+2cj) + ^ sin 2 icos(3M+2w) 


- P 
u 

sin 2 i-^j cosm] 

(B. 31) 

^20 “ ^2 

(^a) [^§ sin 2 isin(2M+2«) 



3 7 

+ sin 2 isinfM+2cj) - sin 2 isin(3M+2cj) 

+ ^ sin 2 i~|| sinM] (B.32) 


For the integrations the approximation M+w = n has been used. 



APPENDIX C 


COMPUTATION OF SECOND ORDER RADIAL ORBIT ERROR 
From equation (4.82) the second order radial orbit error is 


II 

{*4 

a 

u 

-(l-u) 

{ ff Aa 0 dt 

+ a j 

\ If Aa ° dt 


+ 

(1-u) j 

[ |f ® rAtdt 

- a J 

( If arAtdt 


+ 

(l-u) J 

[ |n M r Atdt 
1 3M 

- a j 

[ |n M r Atdt 
1 3M 


+ 

(l-u) J 

[ I 5 u r Atdt 

1 dcj 

- a j 

| cj r 4tdt 

1 OCJ 

(C.l) 


which can be written as 


Ar C2 = (l-u) Aa C2 - aAu r>2 


(C . 2) 


with 


Aa G2 = I a + I a + I a + I 3 (C . 3) 

b a e M w 

aAu G2 = al u + al u + al u 4 al u (C.4) 

b a e M u 


First, the second order effects Aa G2 are going to be computed. 
Since the main effects in a are the ones due to J 2 , a can be 
approximated by a 20 which is given by (B.27). Then the partials 
contained in the integrals of (C.3) are computed to be 
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da 


ZSL 

da 


= ^ = A t + A 2 

3a ’ exp 1 icit da 


with 


A\ - — ^ nJ 2 (^] sin 2 isin(2M+2«) + o(e) 

A 2 = - 6na (- ^ ^ 4tjj 2 |^j sin 2 icos(2M+2«) + o(e) 


where 


3M B 3(nAt) 
da da 


3 n 

2S 4t 


has been used. 


da 


de 


^ ~ 2naJ 2 ^ sin 2 isin(M+2<j) + sin 2 i-~jsi 


inM 


^ sin 2 isin(2M+2cj) j 


d a 2 o _ d a 2 o _ 
dM " d w 


= - 6naJ 2 sin 2 icos(2M+2cj) + o(e) 


Having computed the partials, the integrals of equation (C.3) 
easily computed. If becomes 


la - | (Aj + A 2 )Aa 0 


dt 


- | Aj Aa 0 dt = ^ J 2 Aa 0 sin 2 icos(2M 0 (t)+2cjo(t)) 


(C.5) 

(C.6) 
(C. 7) 

(C.8) 

(C.9) 

(C.10) 
can be 

(C.ll) 

(C.12) 


which can have an amplitude of about 3Aa 0 mm for altimeter satellites 
and can be neglected. 
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- | A 2 Aa 0 dt = 3a (~i ~ Aa 0 At] J 2 sin 2 isin(2M 0 (t)+2cj 0 (t) ) 

(C. 13) 


The integrals I| Ip| Ig are straightforward to compute. They are 
equal to 


I a = 2aJ 2 e r At sin 2 icos(M 0 (t)+2« 0 (t)) + 


(| sin 2 i-|jcosM 0 (t.) - || sin 2 icos(3M 0 (t)+2cj 0 (t.))] 
I a = - 3aJ 2 f^l M r Atsin 2 isin(2M 0 (t)+2cj 0 (t) ) 

I a = - 3aJ 2 [■^ a l 2 r Atsin 2 isin(2M 0 (t)+2cj 0 (t)) 

CJ v 3' 


(C. 14) 


(C . 15) 


(C. 16) 


Then (C.3) becomes 

Aa G2 = - 3aJ 2 ^ Aa 0 + M r + <a r j Atsin 2 isin(2M 0 (t)+2« 0 (t ) ) +Ig 

(C. 17) 

During the integrations (C.12) - (C.16) the following 

approximations have been used 


w+M = « 0 (t) + M 0 (t) (C. 18) 

and 

6+M = n (C. 19) 


The computation of Au G2 requires an expression for u. From 
(4.47) we have 


u = 


du 

dt 


ecosM - eMsinM 


(C.20) 
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A very good approximation to u can be obtained if we approximate e, 
e, M and M of (C.18) by the corresponding effects due to J 2 . Then 

e s e 20 


e - e Q + e 2 o 


M a M 0 + M 2 o s n + M 20 (C.20) • 

cosM = cosM 0 (t) - M 20 sinM o (t) 
sinM s sinM 0 (t) + M 20 cosM 0 (t) 

where M 20 represents only the periodic perturbations due to J 2 such 
that M 0 (t) = M-M 20 - Similarly M 20 contains only the periodic rates 
due to J 2 . Retaining the osculating value of e whenever it multiplies 
M 20 we obtain 


u = — e 0 nsinM 0 (t) + e 20 cosM o (t) - e 20 nsinM 0 (t) 


- eM 20 sinM 0 ft) - enM 20 cosM 0 (t) (C.21) 

In (C.21) all the terms containing products of J 2 perturbations and/or 
their rates have been neglected. The values of e 20 , e 20 , M 20 and M 20 
can be found in Appendix B. Using these values and carrying out all 
the necessary operations we obtain the following form 


u s — e o nsinM 0 (t) - ^ nJ 2 sin 2 isin(M+M 0 (t.)+2co) (C.22) 


Then, consistent with the approximation (C.18) we have 


u 3 - e 0 nsinM 0 (t) - | nJ 2 |^j sin 2 isin(2M 0 (t)+2w 0 ( t) ) (C.22) 


We can see that (C.23) does not contain any singular elements. 
As a matter of fact it is independent of the osculating eccentricity 
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which was cancelled out in the product eM 20 . Using (C.23) the 
partials contained in the integrals of (C.4) become 



+ U 2 


+ U 3 


(C. 24) 


where 


Ui = | 2 eo sinM 0 (t) 


(C . 25) 


U 2 


3 n 2 
2 a 6 


AtcosM 0 (t) 


(C.26) 


U 3 = - 5 nJ 2 (^) 2 (-| J At|sin z icos(2M 0 (t)+2cjo(t) ) 


(C . 27) 


di\ 

9e 


0 


(C.28) 


= - 5 nJ 2 |^j sin 2 icos(2M 0 (t)+2cj 0 (t) ) 


(C . 29) 


9 ii 
do? 


2 

- 5 nJ 2 sin 2 icos(2H 0 ( 1)^2 cj 0 ( l) ) 


(0 . 30 } 


Having computed the partials, the computation of the integrals in (C.4) 
becomes straightforward. 1^ is equal to 


1“ = -a f (Uj + U 2 + U 3 ) Aa 0 dt 


(C.31) 


where 


-a J Ui Aa 0 dt = 0 + 0(e) (C.32) 

- a | U 2 Aa 0 dt = ae 0 [~| 2Aa 0 At]sinM 0 (t) 


(C.33) 
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-a | U 3 Aa 0 dt = | aJ 2 (^] [-| ^ Aa o At]sin 2 isin(2M 0 (t)+2« 0 (t)) 

(C.34) 

The integrals alJJ and al“ are 

al„ = - | aJ 2 (^) VAtsin 2 isin(2M 0 (t)+2« 0 (t)) (C.35) 

al^ = - ^ aJ 2 [%*] cj r Atsin 2 isin(2M 0 (t)+2« 0 (t)) (C.36) 

Then aAu G2 becomes 
3 

aAu G2 = - g e 0 nAa 0 AtsinM 0 (t) 

- ^ aJ 2 [^] ~ Aa 0 +M r +w r j Atsin 2 isin(2M 0 (t)+2cj 0 ( t) ) 

(C.37) 

The final expression for the second order radial orbit error given by 
(C.2), with the term 1-u ~ 1, is 

3 

Ar G2 = ^ e 0 nAa 0 AtsinM 0 (t) 

- ^ aJ 2 ^ Aa 0 + M r +cj r j Atsin 2 isin(2M 0 (t)+2w 0 (t) ) 

+ 2aJ 2 e r At sin 2 icos(M 0 (t) t-2w 0 (t) ) 

+ [| sin 2 i-|] cosM 0 (t) - || sin 2 icos(3M 0 (t)+2w 0 (t))| 

+ J 2 (^) 2 Aa 0 sin 2 icos(2M 0 (t)+2 Mo (t)) (C.38) 

All the above expressions are computed on the reference orbit. 



APPENDIX D 


MULTIPLE ANGLE TO SINGLE ANGLE TRIGONOMETRIC TRANSFORMATION 

The cosine and sine of an angle kv-mu can be transformed to 
functions of cosines and sines of single angles v and u as follows 

cos(kv-mu) = coskvcosmu + sinkvsinmu (D.l) 


sin(kv--mu) = sinkvcosmu - coskvsinmu (D.2) 

From trigonometry we generally have 


cosnx = Re £ f n ]j s cos n s x sin s x 
s-o Isl 


= I ( n ) (-1) s / 2 cos n s x sin 5 x ''P.3) 

S-0 Isl 


where s is even and 


fn if n even 
ln-1 if n odd 


Similarly 


sinnx = Re £ [ n ]j B l cos n s xsin s x 

s =o 

= X [ n ] C — 1 ) 2 cos ,, “ R x sin s x (D.4) 

where s is odd and 
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. fn if n 
n ” in-1 if n 


odd 

even 


Then the four trigonometric products of (D.l) and (D.2) are 


k' m' fl,} f" tti ^ s + t 

coskvcosrau = s I o t I 0 [ t J(-l) 2 F u> 


(D.5) 


where s, t, k' , m' are even 


sinkvsinmu - -J, J, (J) (") 


(D.6) 


where s, t, k' , m' , are odd 


sinkvsinmu = - J i JJ*) (j](-l) “2 ‘ F ux 


CD. 7) 


where s, k' are odd and t, m' are even. 


coskvsinmu = I I P*] [™] (-1) — - ^ — 1 F uv 

g-0 


CD. 8) 


where s, k' are even and t, m' are odd. In all the equations 
(D.5)-(D.8) F uv is 


F uv = cos k s v sin n v cos m f u sin’^u (D.9) 


Equation (D.5)-(D.8) are valid when both k and m are positive 
integers. When both are negative (D.5) and (D.6) remain the same 
while (D.7) and (D.8) change. When one integer is negative then (D.6) 
and (D.7) or (D.8) change. For k negative (D.6) and (D.7) become 


sinkvsinmu = (t) (_1)2 ^ F “' 


sinkvcosmu = f J i JJ'g') (“) (_i)=±i±i p u> 


(D. 10) 


2 


(D.ll) 
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where now F uv contains the absolute value of k. Then 
(D.6), (D.10) and (D.7), (D.ll) we can write for any k 


sinkvsinrau - ^ ' t | 0 ( '£') (*) (-l)V, 


with p = 


f3s+t 


| s+2 
2 


k>0 


k<0 


with c = 


3s+t+l 


s+t+1 


1 f ,kl ] 

Mr 

t-.-o l s J 

ItJ '• 

100 


k<0 



: uv 


Then equation (D.l) is written 
1 k I 


• k i m r | v n ( m ) a , , 

cos(kv-mu) = Z I , (~l) p cos ,kl 3 v sin s v cos m ' 

S-0 t-oK s 


where s, t have the same parity. 
Similarly (D.2) is written 


sin(kv-mu) = X I f™| (-l) c cos* k * s v sin s v cos m 1 

S — O t -- O ^ S ' \ L * 


where s, t have different parities. Equations (D.16) and 
valid for positive in and any k. 


combining 
(D. 12) 

(D. 13) 

(D. 14) 
(D. 15) 

i sin t u 

rn.ifii 

u sint-u 
(D. 17) 

(D.17) are 


APPENDIX E 


DETERMINATION OF CROSSOVER INTERSECTIONS 

As described in Chapter 8 crossover intersections occur at times 
tj and tj that correspond to particular points of the ground track 
with geocentric latitude 4> and longitude X. The determination of 
crossover times and geocentric coordinates can in principle be made 
sequentially by using the satellite ephemeris. The satellite earth 
fixed position at each epoch has to be compared with all the positions 
at following epochs so that to locate the vicinity of a crossover point. 
This process though is extremely time consuming and requires a lot of 
computer memory. 

A more efficient determination can be made using an analytical 
approach. In such a method the secular perturbations of the satellite 
are used to determine approximate times and the geocentric latitude 
and longitude for each intersection. Then, using these approximate 
estimates the exact intersection times and coordinates can be found 
by simply using the precise ephemeris and interpolating. For the 
analytical determination the times and longitudes of all the equator 
crossings for each arc have to be known. Very accurate equator 
crossings can be easily obtained using the precise ephemeris. An 
approximation to the equator crossings can be obtained by simply 
considering that the argument of the satellite, when it crosses the 
equator, takes the value 


(« + f)’ = [2(n - 1) + i ] n 


(E.l) 


where n is the revolution number and i indicates whether the satellite 
crosses the equator towards the north hemisphere (i = 0) or towards 
the south hemisphere (i = 1). The times of equator crossings can be 
found from 


(« + *)A = « 0 + fo + (* + f)tj (E.2) 

which can be combined with (E.l) to give 
t i = [2(n-l)+i]7T _ op + fp 

"> E A 4- -P A 4. f 
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where «, f are the secular rates of the corresponding elements. In 
order to find the longitudes of equator crossings, we consider the 
following 


X = (a - Q) + (0 - 0) 

. , ^ _ sin(w+f)cosi 

sin(a - Q) = 1 r 

cos<t» 

Then at the equator sin(a-Q) = 0 and 


(E.4) 

(E.5) 


(a - Q),!, = s [2(n-l)+i]rr (E.6) 

where s is equal to 1 for posigrade orbits and -1 for retrograde 
orbits. Using (E.4) and (E.6) we obtain 


K = (0 o - 0 o) + (6 - 9 < e + st2(n-l)+i]rr (E.7) 

Or, if we substitute for the time from (E,3) 

X„ = - e 0 ) + [2(n-l)+i]n(s + ^|] - («„ + f„ (E.8) 

v <j+f ' «+ f 


When the inclination is greater than yu” (retrograde orbit) the 
longitude separation between two successive crossings is greater than 
180”. Therefore there is always one intersection between any two 
north and any two south half revolutions and possibly there are two. 
In order to have a second intersection between any two half 
revolutions, k and j, the following condition needs to apply. The 
modulus (with respect to 180”) of the longitude separation between 
the k and j ascending equator crossings has to be smaller than the 
modulus of the longitude separation between the ascending and 
descending equator crossings of the k revolution. Using (E.8) the 
first modulus is 


mod (AX j 
for j > k. 


k , 7 t ) = mod(2( j-k)Tr tt) 

CJ+f 

The second modulus becomes 


0-0 

a+i’ 


Q-8 

W+f 


(E.9) 


modCAX^ 1 , 


7T ) = mod (rr 


(E. 10) 
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Then if (E.9) is smaller than (E.10) there are two intersections 
between the k and j half (north or south) revolutions. 

To find the longitude of a crossover intersection between any k 
and j revolutions in the north hemisphere the following relationship is 
used. 


x jk = (K + X j + K + x j)/4 (E.ll) 

if 

. 0 1 O 1 

| Ak * Xk X j | < 7T (E. 12) 

When the expression in (E.12) is greater than tt then Xj k needs to be 
increased by n. The equator crossing longitudes in (E.ll), (E.12) 
need to satisfy 


> « x k 


(E.13) 


sXj > sXl 


If this is not the case, 2-n has to be added to the appropriate 
longitudes to satisfy (E.13). In the southern hemisphere the 
longitude of a crossover intersection between any k and j revolutions 
can be obtained by 


x jk ~ ( x k + x ] + x k+i + x j+i)/ 4 


(E . 14) 


with conditions similar to the ones for (E.ll) being valid. When more 
than one intersection (north or south) is formed between two 
revolutions, the longitude of the second intersection is simply offset 
by 7i from the longitude of the first one. 

In order to find the time of the intersection we use 


tan(w + f) 


tanA 


cosi 


(E. 15) 
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where A is the longitude separation between the intersection and the 
equator crossing immediately preceding the intersection. This 
longitude separation for the k revolution is 

a! = A - s(Q - 0)r’ (E. 16) 

where 

A= |x Jk -Aj| (E. 17) 


is the corresponding longitude separation at time tj. jE . The second 
term of (E.16) is the additional longitude increment corresponding to 
the time r{. that is required for the satellite to move from the equator 
to the intersection point. Considering the above, (E.15) becomes 


tant(£, + f)r£] = tan[t - »jn-9)r} ] 
k cosi 


(E . IS) 


Equation (E.18) can be solved iteratively for until convergence. 
The same determination can be made for -rj. Then the crossing times 
tj. and tj are 


tj = t! 


k,E 


(E. 19) 


t; - t;, £ + t 


(F..20) 


When a second crossing between two half revolutions is formed the 
corresponding times are 


* i _ 4. i , i 
^k ^k , E + T j 


+ > _ + 1 , 1 
t j " tj.E + T k 


(E.21) 

(E.22) 


The latitude of the crossing is simply found by 


♦ jk = sin 1 [sin(w 0 +f 0 +(«+f)t j)] 


(E . 23) 
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The same value is obtained by using tj instead of tj. The total 
number of crossings that are formed from N complete revolutions is 


M = N(N-l) + (1-s) l J 6 ik 

j=i k=i J 


(E.24) 


where 


6 


jk 


1 if (E. 9) < (E . 10) 
0 otherwise 


(E . 251 


The analytic approach for determining the crossing locations is 
not very accurate for two reasons. First, the full perturbations of 
cj,f,fl have to be used and not only their secular counterparts. 
Second, the assumption of a circular orbit, which is not valid, has 
been used. So the crossover solutions need to be improved by 
computing corrections to the analytic solutions. These corrections can 
be obtained by a nonlinear interpolation scheme of the precise times 
and coordinates taken from the precise ephemeris before and after 
the crossing. 



